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ABSTRACT 


Time-dependent  numerical  models  are  developed  to  predict  the  temporal  and 
cross-shore  variations  of  the  free  surface  elevation  and  the  cross-shore  and  alongshore 
fluid  velocities  in  the  swash  and  surf  zones  under  obliquely  incident  waves.  The  as¬ 
sumption  of  shallow  water  with  small  incident  angles  and  slow  alongshore  variations 
are  made  to  reduce  computational  efforts  considerably  and  to  eliminate  difficulties 
associated  with  lateral  boundary  conditions.  These  assumptions  enable  the  mod¬ 
els  to  compute  the  cross-shore  fluid  motion  separately  from  the  alongshore  motion. 
The  numerical  models  allow  gradual  alongshore  variations  of  the  bathymetry  and 
the  incident  regular  or  irregular  waves  specified  at  the  seaward  boundary. 

Two  numerical  models  are  developed  in  this  study.  The  first  model  is  a  two- 
dimensional  (2D)  model  that  neglects  the  vertical  variations  of  the  cross-shore  and 
alongshore  velocities.  As  a  result,  this  model  neglects  the  dispersion  due  to  the 
vertical  variations  of  the  horizontal  velocities  and  predicts  only  the  depth-averaged 
cross-shore  and  alongshore  velocities.  The  second  model  is  a  quasi  three-dimensional 
(3D)  model  that  assumes  a  cubic  profile  for  the  horizontal  velocities  and  includes 
the  dispersion  terms  due  to  the  vertical  variations  of  the  horizontal  velocities.  Two 
additional  equations  for  cross-shore  and  alongshore  momentum  flux  corrections  are 
derived. 

These  numerical  models  are  compared  with  available  laboratory  and  field 
data  for  planar  beaches.  To  assess  the  importance  of  the  dispersion  terms,  both 
models  are  compared  with  the  same  data.  Both  models  are  shown  to  be  capable 
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of  predicting  the  cross-shore  variations  of  wave  height,  setup  and  runup  for  regu¬ 
lar  waves  and  root-mean-square  wave  height  for  irregular  waves.  The  dispersion 
effects  on  wave  height  and  setup  are  shown  to  be  minor.  The  models  predict  large 
cross-shore  and  alongshore  velocities  near  the  shoreline  for  irregular  waves.  The 
causes  of  these  large  velocities  are  examined  but  can  not  be  ascertained  for  lack  of 
velocity  data  near  the  shoreline.  For  regular  waves,  the  2D  model  with  the  bottom 
friction  factor  calibrated  previously  for  swash  oscillations  predicts  the  magnitude 
of  longshore  current  but  cannot  reproduce  the  longshore  current  profile.  For  ir¬ 
regular  waves,  the  longshore  current  is  predicted  fairly  well  by  both  models.  The 
3D  model  improves  the  prediction  of  the  longshore  current  profile  for  regular  waves 
significantly.  This  implies  that  the  dispersion  effects  on  the  longshore  current  are 
significant  for  regular  waves  but  secondary  for  irregular  waves.  The  3D  model  is 
also  shown  to  predict  the  vertical  shape  of  the  longshore  current  inside  the  surf  zone 
but  not  outside  the  surf  zone.  The  vertical  variation  of  the  longshore  current  for 
irregular  waves  is  shown  to  be  small. 

The  3D  model  is  also  compared  with  field  data  for  a  barred  beach.  The 
model  underpredicts  the  root-mean-square  wave  height  in  the  bar  trough  region. 
The  field  measurements  of  longshore  currents  generally  indicate  a  broad  peak  in  the 
bar  trough  region.  Under  the  assumption  of  alongshore  uniformity,  the  model  cannot 
explain  these  broad  peaks.  The  small  alongshore  variation  of  wave  setup  induced 
by  a  small  alongshore  variation  of  obliquely  incident  irregular  waves  is  shown  to 
significantly  modify  the  driving  force  and  longshore  current  profile  in  the  bar  trough 
region  where  the  cross-shore  gradient  of  the  alongshore  radiation  stress  is  small. 
As  a  result,  the  longshore  current  profile  in  the  bar  trough  region  is  sensitive  to 
the  alongshore  variability.  On  the  other  hand,  for  planar  beaches,  the  longshore 
current  profile  is  shown  to  be  insensitive  to  the  alongshore  variation  of  obliquely 
incident  waves.  This  may  explain  why  existing  longshore  current  models  based  on 
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the  assumption  of  alongshore  uniformity  were  regarded  to  be  adequate  before  their 
comparisons  with  the  barred  beach  data. 


Chapter  1 


INTRODUCTION 


Obliquely  incident  waves  breaking  on  a  beach  generate  longshore  current 
that  in  turn  drives  longshore  sediment  transport.  Knowledge  of  longshore  sediment 
transport  is  important  for  the  design  of  shoreline  erosion  mitigation  measures  such  as 
sand  bypassing  and  beach  nourishment.  Longshore  sediment  transport  mostly  takes 
place  within  the  region  where  the  waves  deform  after  breaking  and  the  region  where 
a  thin  sheet  of  water  rushes  up  and  down  the  beach,  known  as  the  surf  zone  and 
swash  zone,  respectively.  Therefore,  a  quantitative  understanding  of  surf  and  swash 
hydrodynamics  on  beaches  under  obliquely  incident  waves  is  required  for  modeling 
of  sediment  transport  in  these  regions. 

Field  and  laboratory  measurements  on  the  distribution  of  longshore  sediment 
transport  across  the  surf  zone  indicated  that  the  distribution  was  generally  bimodal 
with  peaks  in  the  swash  and  breaker  zones  {e.g.^  Bodge  and  Dean  1987;  Kamphuis 
1991).  Bodge  and  Dean  (1987)  observed  that  the  relative  significance  of  the  peaks 
shifted  from  the  breaker  zone  peak  to  the  swash  zone  peak  as  the  surf  varied  from 
spilling  to  collapsing  conditions.  They  found  that  longshore  sediment  transport  in 
the  swash  zone  might  account  for  at  least  5%  to  over  60%  of  the  total  drift. 

There  are  two  different  approaches  to  study  the  surf  zone  hydrodynamics, 
namely,  time-averaged  and  time  dependent  models.  Time-averaged  models  split  the 
fluid  motion  into  a  ‘wave’  part  and  a  ‘current’  part  where  it  is  assumed  that  the 
wave  motion  is  already  known.  On  the  other  hand,  time-dependent  models  account 


for  wave  propagation  and  there  is  no  need  to  differentiate  between  the  wave  and 
current  parts  of  the  motion.  This  is  significant  when  dealing  with  the  complicated 
wave  motion  in  the  surf  and  swash  zones,  although  this  type  of  model  requires  more 
computational  effort. 

Since  the  introduction  of  the  concept  of  radiation  stresses  by  Longuet-Higgins 
and  Stewart  (1960),  time-averaged  models  have  been  applied  widely  to  predict  the 
cross-shore  variations  of  the  wave  height,  setup,  longshore  current  and  sediment 
transport  rate  in  the  surf  and  swash  zones  along  a  straight  shoreline  {e.g.,  Battjes 
1988).  These  models  are  relatively  simple  and  useful  for  practical  applications  except 
in  the  swash  zone.  Contrary  to  the  field  and  laboratory  data,  the  models  based 
on  the  time-averaged  momentum  and  energy  equations  predict  that  the  longshore 
current  and  sediment  transport  rate  approach  zero  at  the  mean  water  shoreline. 
These  models  also  underpredict  the  wave  height  and  setup  in  the  swash  zone  (Cox 
et  al.  1994).  Furthermore,  swash  measurements  on  natural  beaches  (e.^.,  Guza  and 
Thornton  1982;  Holman  and  Sallenger  1985)  show  the  dominance  of  low-frequency 
wave  motions  unless  incident  wind  waves  collapse  or  surge  on  foreshore  slopes. 

The  main  aim  of  this  study  is  to  develop  time-dependent  models  in  an  attempt 
to  predict  the  time-dependent  and  time-averaged  hydrodynamics  in  the  surf  and 
swash  zones  as  mentioned  above. 

1.1  A  Brief  Overview  of  Related  Works 

To  explain  a  large  volume  of  sediment  transport  in  the  swash  zone,  Thornton 
and  Abdelrahman  (1991)  proposed  an  analytical  model  by  including  a  specified 
standing  wave  of  low-frequency  interacting  with  an  obliquely  incident  wave  of  a  wind- 
wave  frequency.  As  a  result  of  the  interaction  between  incident  and  low-frequency 
waves,  a  non- vanishing  longshore  current  and  sediment  transport  rate  were  found  at 
the  shoreline.  This  model  indicates  the  importance  of  low-frequency  standing  waves 
in  the  swash  zone  that  need  to  be  estimated  or  predicted  separately.  To  predict  the 
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interaction  between  incident  and  low  frequency  waves  in  a  more  realistic  manner, 
a  time-dependent  numerical  model  will  be  required  although  it  will  require  more 
computational  effort. 

Ryrie  (1983)  developed  a  time-dependent  numerical  model  for  longshore  fluid 
motion  generated  by  obliquely  incident  monochromatic  waves  with  a  small  angle  of 
incidence  along  a  straight  shoreline  with  a  plane  slope.  The  model  is  based  on  the 
two-dimensional  shallow  water  equations.  The  cross-shore  and  alongshore  motions 
are  decoupled  using  a  new  time  variable,  referred  to  as  pseudo-time,  which  is  a 
function  of  the  phase  speed  and  incident  angle  of  the  incident  waves  at  the  break 
point.  The  numerical  model  was  not  compared  with  any  data.  Recently,  Brocchini 
and  Peregrine  (1996)  examined  the  integral  flow  properties  of  the  swash  zone  using 
the  model  of  Ryrie  (1983)  and  the  standing  wave  solution  of  Carrier  and  Greenspan 
(1958). 

Kobayashi  et  al.  (1987,  1989)  and  Kobayashi  and  Wurjanto  (1992)  developed 
and  evaluated  a  time-dependent,  one-dimensional  numerical  model  for  normally 
incident  waves  using  laboratory  and  field  data  on  wave  runup  on  steep  and  gentle 
slopes.  This  numerical  model  was  shown  by  Cox  et  al.  (1994)  to  agree  fairly  well 
(within  about  20%  errors)  with  the  free  surface  measurements  in  the  swash  and 
inner  surf  zones  using  ten  wave  gages,  partially  buried  in  the  sand,  in  a  large  wave 
flume.  Similar  results  were  also  observed  from  the  measurements  by  Raubenheimer 
et  al.  (1995)  using  a  vertical  stack  of  five  runup  wires  supported  parallel  to  a  gently 
sloping  natural  beach.  Karjadi  and  Kobayashi  (1994)  and  Kobayashi  and  Karjadi 
(1994a)  modified  this  model  to  predict  solitary  wave  breaking,  runup  and  reflection. 
The  existing  one-dimensional  model  is  expanded  herein  to  predict  the  cross-shore 
and  temporal  variations  of  the  free  surface  elevation,  the  alongshore  and  cross-shore 
fluid  velocities  in  the  swash  and  surf  zones  under  obliquely  incident  waves.  The 
expanded  numerical  model  is  formulated  unlike  that  of  Ryrie  (1983)  such  that  it  is 
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applicable  to  beaches  of  arbitrary  geometry  under  obliquely  incident  random  waves. 

To  reduce  computational  efforts  considerably  and  to  eliminate  difficulties  as¬ 
sociated  with  lateral  boundary  conditions  for  general  two-dimensional  computations, 
the  assumptions  of  small  incident  angles  and  slow  alongshore  variations  are  made. 
These  assumptions  enable  the  model  to  compute  the  alongshore  velocity  separately 
after  the  free  surface  elevation  and  the  cross-shore  velocity  are  computed  using  the 
one-dimensional  model  along  at  least  three  lines  normal  to  the  shoreline.  On  the 
other  hand,  because  of  these  assumptions,  the  numerical  model  expanded  herein 
cannot  be  used  to  predict  edge  waves  (e.^.,  Bowen  and  Guza  1978;  Huntley  et  al. 
1981),  which  appear  to  play  an  important  role  in  the  generation  of  rhythmic  beach 
morphology  (e.^.,  Holman  and  Bowen  1982).  Furthermore,  rip  currents  (e.y.,  Dal- 
rymple  1978)  and  shear  waves  (Bowen  and  Holman  1989;  Oltman-Shay  et  al.  1989) 
are  beyond  the  capabilities  of  the  expanded  numerical  model. 

1.2  Outline  of  Report 

In  Chapter  2,  the  governing  equations  employed  in  the  expanded  numerical 
model  are  derived  from  the  three-dimensional  continuity  and  Reynolds  equations 
(e.^.,  Rodi  1980)  in  a  manner  similar  to  the  derivation  of  Kobayashi  and  Wurjanto 
(1992)  of  the  approximate  one- dimensional  equations  from  the  two-dimensional  con¬ 
tinuity  and  Reynolds  equations. 

The  approximate  three-dimensional  equations  under  the  assumptions  of  shal¬ 
low  water  waves  on  a  relatively  gentle  impermeable  slope  with  small  angles  of  in¬ 
cidence  are  derived  from  the  three-dimensional  continuity  and  Reynolds  equations. 
The  approximate  three-dimensional  equations  are  then  integrated  over  the  water 
depth  to  obtain  the  depth-integrated  continuity  and  horizontal  momentum  equa¬ 
tions.  In  this  study,  two  numerical  models  are  developed.  The  first  model  is  based 
on  the  approximate  two-dimensional  equations  that  neglect  the  dispersion  terms  due 
to  the  vertical  variations  of  the  horizontal  velocities.  This  model  is  simply  referred 
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to  as  the  2D  model.  The  second  model,  based  on  a  cubic  profile  of  the  horizontal  ve¬ 
locities,  includes  the  dispersion  terms  due  to  the  vertical  variations  of  the  horizontal 
velocities.  For  brevity,  this  quasi-3D  model  is  referred  to  as  the  3D  model. 

The  numerical  procedures  to  solve  the  governing  equations  with  appropriate 
initial  and  boundary  conditions  are  explained  in  Chapter  3.  For  the  2D  model,  the 
governing  equations  for  the  cross-shore  wave  motion  for  small  angles  of  incidence, 
which  are  the  same  as  those  of  the  one-dimensional  model,  are  solved  using  the 
Lax-Wendroff  method  in  the  same  way  as  the  one- dimensional  model  developed  by 
Kobayashi  et  al.  (1987,  1989).  For  the  longshore  motion  in  the  2D  model  and 
for  both  the  cross-shore  and  alongshore  motions  in  the  3D  model,  the  MacCormack 
(MacCormack  1969)  finite  difference  method  is  selected  because  of  its  simplicity  and 
success  in  the  computation  of  unsteady  open  channel  flows  with  hydraulic  jumps 
(Chaudhry  1993). 

In  Chapter  4,  the  2D  numerical  model  is  compared  with  laboratory  data 
of  Visser  (1991)  and  field  experiment  conducted  by  Thornton  and  Guza  (1986)  at 
Leadbetter  beach  in  1980.  The  comparisons  indicate  that  the  2D  model  predicts 
the  maximum  setup  and  runup  fairly  well  although  the  model  does  not  predict  the 
transition  zone  of  constant  wave  set-down.  Also  for  regular  waves  the  model  cannot 
predict  the  longshore  current  profile  since  the  model  does  not  include  the  disper¬ 
sion  due  to  the  vertical  variation  of  the  horizontal  velocities.  For  irregular  waves, 
the  model  predicts  the  root-mean-square  wave  height  and  longshore  current  fairly 
well.  This  implies  that  the  effects  of  the  dispersion  on  the  longshore  current  under 
irregular  waves  are  negligible  as  already  concluded  by  Thornton  and  Guza  (1986). 
The  model  predicts  large  cross-shore  and  alongshore  velocities  near  the  shoreline 
and  very  low-frequency  oscillations  of  the  depth-averaged  alongshore  velocity.  The 
causes  of  these  results  are  also  examined  in  this  chapter.  The  summaries  of  these 
comparisons  were  published  by  Kobayashi  and  Karjadi  (1994b,  1996). 


In  Chapter  5,  the  3D  model  is  compared  with  the  same  regular  and  irregular 
wave  data  as  the  2D  model  to  assess  the  importance  of  the  added  dispersion  terms. 
These  limited  comparisons  indicate  that  the  significant  improvement  due  to  the 
addition  of  the  momentum  flux  corrections  associated  with  the  vertical  variations 
of  the  horizontal  velocities  is  essentially  limited  to  the  cross-shore  profile  of  the 
longshore  current  induced  by  regular  breaking  waves.  The  effects  of  the  dispersion 
terms  on  wave  height,  setup  and  runup  are  shown  to  be  secondary.  The  3D  model 
is  also  shown  to  predict  the  vertical  shape  of  the  longshore  current  inside  the  surf 
zone  but  not  outside  the  surf  zone.  The  summary  of  Chapter  5  will  be  published  by 
Kobayashi,  Karjadi  and  Johnson  (1997). 

Chapter  6  discusses  the  application  of  the  3D  model  to  a  barred  beach.  The 
model  is  compared  with  the  DELILAH  field  data  of  Smith  et  al  (1993).  The  mea¬ 
surements  of  longshore  current  on  the  barred  beach  during  the  DELILAH  experiment 
generally  indicated  a  broad  peak  in  the  bar  trough  region.  The  3D  model  cannot 
explain  this  data  under  the  assumption  of  alongshore  uniformity.  Low-frequency 
components  and  alongshore  variations  of  incident  irregular  waves  are  examined  to 
explain  the  broad  peak  of  the  longshore  current.  The  longshore  current  profile  on 
the  barred  beach  is  then  shown  to  be  sensitive  to  the  alongshore  variability  unlike 
the  longshore  current  profile  on  the  plane  beaches  examined  in  Chapters  4  and  5. 
The  summary  of  Chapter  6  will  be  published  by  Karjadi  and  Kobayashi  (1996). 

The  summary  and  conclusions  of  this  study  is  given  in  Chapter  7. 
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Chapter  2 

MATHEMATICAL  FORMULATION 

The  equations  that  govern  the  obliquely  incident  wave  motions  on  an  imper¬ 
meable  slope  are  derived  and  discussed  in  this  chapter.  The  approximate  governing 
equations  under  the  assumptions  of  shallow  water  waves  with  small  angles  of  in¬ 
cidence  are  derived  from  the  three-dimensional  continuity  and  Reynolds  equations 
(e.g.,  Rodi  1980).  This  derivation  is  similar  to  that  of  the  one-dimensional  equa¬ 
tions  from  the  two-dimensional  continuity  and  momentum  equations  proposed  by 
Kobayashi  and  Wurjanto  (1992). 

The  notations  used  in  the  following  analysis  are  depicted  in  Figure  2.1.  The 
angle  of  wave  incidence  is  9c  and  will  be  assumed  small.  The  other  symbols  shown  in 
this  figure  will  be  explained  in  the  subsequent  derivation.  The  model  is  formulated 
in  such  a  way  that  it  will  be  applicable  for  beaches  with  arbitrary  geometry  under 
obliquely  incident  regular  or  random  waves. 

2.1  Approximate  Three-Dimensional  Equations 

The  three-dimensional  continuity  and  Reynolds  equations  are  (e.g.,  Rodi 


where  the  prime  indicates  the  physical  variables  and  the  summation  convention  is 
used  with  respect  to  repeated  indexes.  The  symbols  involved  in  (2.1)  and  (2.2)  are 
depicted  in  Figure  2.1  and  explained  in  the  following: 


Figure  2.1:  Definition  sketch  for  obliquely  incident  wave  motion  on  an  imperme¬ 
able  slope. 


i'  =  time 

x\  =  horizontal  coordinate  normal  to  the  overall  shoreline  orientation, 
taken  to  be  positive  landward 
x'2  —  horizontal  coordinate  normal  to  x\ 

Xg  =  vertical  coordinate  taken  to  be  positive  upward  with  Xg  =  0 
at  the  still  water  level  (SWL) 
u\  =  horizontal  velocity  in  x'-^  direction 
Ug  =  horizontal  velocity  in  Xg  direction 
Ug  =  vertical  velocity 

p  =  fluid  density,  which  is  assumed  constant 
p'  =  pressure 

g  =  gravitational  acceleration 
Si3  =  Kronecker  delta 

T--  =  sum  of  turbulent  and  viscous  stresses 
u 

Assuming  that  the  viscous  stresses  are  negligible,  may  be  expressed  as 
{e.g.,  Rodi  1980) 

^lj  =  P 

where 

i/(  =  turbulent  eddy  viscosity 

k'  =  turbulent  kinetic  energy  per  unit  mass 

To  simplify  (2.1)  and  (2.2)  together  with  (2.3)  in  shallow  water,  the  dimen¬ 
sional  variables  may  be  normalized  as 


t' 

^  'J't  ’ 

x\ 

(2.4) 

X 

^  T'^gWyi^  ’ 

Xo- 

’  ^  H' 

“i 

Uo 

II 

(2.5) 

(2.6) 

{gH>y/^  ’ 

eygH'yi^  ’ 

h- 

T'igH'yf^ 

fft2 1'J'l  ’ 

H'{gH'y/yT' 

’  H' 

'du\ 


+ 


du'P 

dx\ 


-  Ik'Sij 


(2.3) 
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where  T\  H'  and  Oc  are  the  characteristic  wave  period,  height  and  incident  angle  in 
radian  used  for  the  normalization.  The  parameter  a  defined  in  (2.6)  is  the  ratio  of 
the  cross-shore  length  scale  to  the  vertical  length  scale.  The  normalized  variables 
in  (2.4)  and  (2.5)  are  assumed  to  be  on  the  order  of  unity  in  shallow  water.  The 
normalization  of  v[  and  k'  in  (2.6)  is  based  on  the  turbulence  measurements  in  a 
wave  flume  by  Cox  et  al.  (1994),  which  have  indicated  that  ut  and  k  are  on  the 
order  of  unity  or  less  inside  and  immediately  outside  the  surf  zone,  respectively. 

Substituting  (2.4)  and  (2.5)  into  (2.1),  the  normalized  continuity  equation 


can  be  written  as 


dxr  ^  dx2  dx3 


0 


(2.7) 


Substituting  (2.4),  (2.5),  and  (2.6)  into  (2.2)  and  (2.3),  the  horizontal  momentum 
equations  are  given  by 

■  1 


dui  ,  dui  2  ,  ^“1 


dp  d 
dxi  dx^ 


—z  2ui 


dui 

dxi 


2  k 

3  a 


.^lA. 

cr^  dxo 


f dui  ^  du2\ 

d 

\5x2  '  ^^1/. 

dxz 

' du\  1  duz' 
'  dxz  0-2  dxx 


and 


du2  du2  ,  „2  dU2  ,  du2 

dt  dxi  dx2  dxz 


dp  1 


d 


d 


dx2 


ei 


2ut 


du2 

dx2 


dx2 
2k 
3  cr 


+ 


dxi 

d 


dxs 


(  du2  du\ '' 

dx2j 

( du2  1  duz^ 
V5x3  cr^  dx2. 


whereas  the  vertical  momentum  equation  is  given  by 

1  ( duz  duz  ,  „2  duz  ,  duz\  d 

+  wi— +  ^>27; - hw3- 


dt 


1 


dx\ 

d 


dx\ 


Vt 


dx2 
1  duz 


+ 


dxz 

dui\ 


cr^  dxi  dx 


+ 


dxz 
cr^  dx2 
dxz 


{p  +  Xz) 

1  duz 


Vt 


—z  2u, 


(T^  dx2 
duz 


-I- 


dxz 


dxz) 

2  k 

3  <7 


(2.8) 


(2.9) 


(2.10) 


Under  the  assumptions  of  >  1  and  0^  <  1  for  shallow  water  waves  with 
small  angles  of  incidence,  the  terms  on  the  order  of  and  smaller  in  (2.7), 
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(2.8),  (2.9),  and  (2.10)  may  be  neglected.  The  normalized  continuity  equation  (2.7) 
may  thus  be  simplified  as 

=  0  (2.11) 

OXi  UXs 


The  normalized  horizontal  momentum  equations  (2.8)  and  (2.9)  with  i  =  1  and  2 


may  also  be  simplified  as 


dui  dui  dui 


(2.12) 


Ti  =  Vt 


and  the  normalized  vertical  momentum  equation  (2.10)  is  simplified  as 


P  +  X3  +  -  - 
6  (J 


(2.13) 


(2.14) 


The  symbols  r]'  and  d'  in  Figure  2.1  represent  the  free  surface  elevation  above 
SWL  and  the  water  depth  below  SWL,  respectively.  Their  normalized  counterparts 


are  given  by 


(2.15) 


where  R'  is  the  characteristic  wave  height  used  for  the  normalization. 

To  further  simplify  the  governing  equations,  the  approximate  three-dimen¬ 
sional  continuity  and  momentum  equations  derived  above  will  be  integrated  over 
the  water  depth  using  the  kinematic  and  dynamic  boundary  conditions  at  the  free 
surface  and  bottom.  The  kinematic  free  surface  and  bottom  boundary  conditions 
require  that  the  total  derivative  of  any  fixed  or  moving  surface  with  respect  to  time 
is  equal  to  zero  on  the  surface  {c.g.,  Dean  and  Dalrymple  1984).  The  kinematic  free 
surface  boundary  condition  is  expressed  as 

^7^  1  ^7  I  _  n  (2.16) 


+  “'23:7-“3  =  0 


whereas  the  kinematic  bottom  boundary  condition  is 

,  dd'  ,  dd!  , 


at  X3  =  77' 


at  *^3  d 


(2.17) 
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(2.18) 


Substitution  of  (2.4),  (2.5)  and  (2.15)  into  (2.16)  and  (2.17)  yields 


^'<1  /i2  n 


at  xs  =  Tj 


and 


dd  2  , 

Ui- - h  6'^«2-5 - i-U3 

UXi  0X2 


0  at  X3  =  —d 

For  61  ^  I,  these  boundary  conditions  may  be  simplified  as 

dr]  di)  A  + 

^  +  ui^ - U3  =  0  at  X3  =  »? 

ot  Ox  I 

and 


Ui 


dd 

dxi 


-)- 1^3  —  0  at  X3  —  d 


(2.19) 


(2.20) 


(2.21) 


For  the  dynamic  boundary  conditions  at  the  free  surface  X3  =  t),  the  tangen¬ 
tial  and  normal  stresses  at  the  free  surface  are  assumed  to  be  zero,  which  can  be 
shown  to  result  in 

n  =  0  with  2  =  1,2  at  X3  =  rj  (2.22) 


and 

p-)-  —  =  0  atx3  =  77  (2.23) 

Bo¬ 
as  explained  in  the  following. 

Figure  2.2  shows  the  surface  stresses  acting  on  an  infinitesimal  fluid  element 
with  area  ABC=dS,  area  OBC=dS'i,  area  OAC=d52,  and  area  0AB=d53.  The 
total  stress  on  the  surface  dS  is  S'  =  (Si,E2,S3)  and  the  unit  normal  vector 
perpendicular  to  dS  is  n'  =  {n'i,n'2,n'^). 

According  to  Gauss  theorem  (e.p.,  Greenberg  1988): 


A.  n'  dS  for  any  A 

With  A=constant,  the  equation  becomes 


(2.24) 


a 


n'dS  =  0 
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where 


S  =  sum  of  dS, dS\,dS2,  and  dSs 

n'=  positive  outward  from  dS,dSi,dS2,  and  dSz 


Hence 

n'  dS  -  el'  dSi  -  el'  dS2  -  el'  dSs  =  0 
from  which  dSi  =  n[  dS,  dS2  =  dS^  and  dS^  =  n'^  dS. 

Using  these  relationships,  the  force  balance  on  fluid  element  OABC  is  given  by 

Xj  —  direction  :  ~  p')dSi  —  T[2dS2  —  +  Ti\dS  =  0 

x'2  —  direction  :  —r'^idSi  —  —  p')dS2  —  r^^dSs  +  =  0 

Xg  —  direction  :  — Tgjd^i  —  T'^2dS2  —  (^33  “  p)dSz  +  Egd^  =  0 
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which  can  be  rewritten  as 


=  (Tjj  P  )n.j  -|-  T\2^2  "b  ^13^3 

^2  ~  ”^21^1  ~b  (^22  P  )^2  "b  '^23^3 
^3  ~  "^31^1  ”b  ^32^2  "b  {'^33  ~~  P  )^3 


VF' 


with  «'  =  |y^  surface  jP'(xj,  X2,  X3,  t')  =  0 

For  the  free  surface  F'  =  x'^  -  p'{x\,x'2,t')  =  0,  the  gradient  of  the  free  surface, 
VF',  is  given  by 


dx'-^  ’  8x2  ’ 


<7  dx\  ’ 


and 


|vr|  = 


SnV  .S_(Sr,y 

'  n 


+  1 


^  1^(72  \^5xi/  a'^  \dx2^ 

Hence,  under  the  assumptions  of  cr^  >  1  and  <  1,  the  unit  vector  ft'  is  given  by 

1  dp  Oc  dp 


■'  =  (ni,ra2>^3)  =  (- 


-,1) 


a  dx-i  ’  <7  dx2 

The  conditions  of  zero  tangential  and  normal  stresses  at  the  free  surface  thus  lead 
to 


0  = 
0  = 
0  = 


1  dp  ' 
a  dx\^ 
1  dp  ' 
cr  dxx^ 
1  dp 
a  dxi 


{<x-p')  + 


Oc  dp 

£7  dX2 


'^12  "b  "^13 


^21  + 


'^31  + 


Oc  dp 
cr  8x2 
1  dp 
cr  dx' 


^  ('^22  “  P')  +  "^23 
^  '^32  +  ("^33  ”  P') 


Normalization  of  the  above  equations  using  (2.3),  (2.4),  (2.5),  (2.6)  and  (2.15)  yields 


0  = 


1  dp 
£7  dx\ 


-  2vt 


du\ 

dxx 


- k  —  crp 

3 


+ 


Oc  dp\  Oc  ( dux  du2 


a  8x2)  cr  ^  \dx2  '  dxx 


+ 


/  dux  1  du3' 

+  I  - 1 - 

\  8x3  cr^  dxx  ^ 


(2.25) 
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0  = 


I  dr]  \ 
a  dxi  J  a  ^ 


'  du2  dui^ 
dxi  dx2, 


+ 


6c  dy 

(7  dxo 


ei 


2vt 


du2 

dx2 


- k  —  crp 

3 


(du2  1  duz^ 

I  -7  I  2~a 
\  oxz  cr^  0X2 , 


(2.26) 


0  = 


1  drj 
a  dx\ 


Vt 


1  duz  dui 


+ 


cr^  dxi  dx 


+ 


+-  [217^ 


6  c  drj  \ 

a  dx2  J 

dxz. 


6c>^t 


1  duz  du2 


+ 


(T^  dx2  dx 


- k  —  ap 

3 


(2.27) 


Neglecting  terms  of  the  order  {a  ^),  (61)  and  smaller,  (2.25),  (2.26),  and  (2.27)  may 
be  simplified  as 

1  drj 


0  = 
0  = 
0  = 


a  dxi ) 

1  dy  \ 
adx2  J 

-k  —  ap 


2 ,  \  (dui 

dxz. 


--k  -ap]+i7t 


(2.28) 

(2.29) 

(2.30) 


Equation  (2.30)  corresponds  to  (2.23),  while  substitution  of  (2.30)  into  (2.28)  and 
(2.29)  gives  the  normalized  boundary  conditions  at  the  free  surface  xz  —  r]  given  by 

(2.22)  with  (2.13). 

Integration  of  the  vertical  momentum  equation  (2.14)  with  respect  to  xz  using 

(2.23)  yields 

p  =  ri-xz-‘^  (2.31) 


which  shows  that  the  pressure  is  approximately  hydrostatic  below  the  free  surface 
elevation  77  where  k  is  on  the  order  of  unity  or  less  and  a  is  relatively  large  to  satisfy 
1.  Substituting  (2.31)  into  (2.12),  the  horizontal  momentum  equations  are 
rewritten  as 

dui  .  dui  dp  .  dTi  ^2  32) 


dui 


dui 

IT  +  3X3 


dp  ^ 
dxi  dxz 


with  z  =  1, 2. 
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The  approximate  continuity  equation  (2.11)  and  the  approximate  momentum 
equations  (2.32)  are  integrated  vertically  to  further  simplify  the  governing  equations 
in  the  following. 


2.2  Depth-Integrated  Equations 

To  reduce  computational  efforts  significantly,  the  approximate  three-dimen¬ 
sional  continuity  and  momentum  equations  are  integrated  over  the  water  depth.  The 
depth-integrated  continuity  equation  is  obtained  by  integrating  (2.11)  with  respect 
to  X3  from  the  bottom,  which  is  assumed  to  be  fixed  and  impermeable,  to  the  free 
surface.  Applying  Leibniz  rule  {t.g.,  Greenberg  1988)  and  the  boundary  conditions 
given  by  (2.20)  and  (2.21),  the  depth-integrated  continuity  equation  is 


d  n 
dx\  J-d 


rv 

/  Ui  dxs  — 
-d 

dr] 

— Ul  -  U3 
uXi 

— 

v. 

xz=r}^ 

d{-d) 

dxi 


-l-  U3 


=  0 


J  xz=—d 


=-i  (2-21) 


=0  (2.20) 


and  can  be  rewritten  as 

dh  dU^h  ^ 

with 

h{t,  Xi.Xi)  -  7/(t,  Xi,  X2)  -f  d{Xi,X2) 


(2.33) 


and 

U^h  =  f^uidxs  (2.34) 

The  depth-integrated  horizontal  momentum  equations  are  derived  in  the  fol¬ 
lowing.  By  use  of  (2.11),  the  horizontal  momentum  equations  (2.32)  can  be  rewritten 


as 

du\ 

duj 

duiUs 

dr]  ^ 

~di 

+  ai. 

H - =  - 

dxs 

dxi  dxz 

and 

du2 

du\U2 

du2U2 

dr]  djj^ 

dxy 

dxs 

dx2  dxz 
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Integrating  the  above  equations  with  respect  to  X3  from  the  bottom  to  the  free 
surface  and  applying  Leibniz  rule  and  the  boundary  conditions  given  by  (2.20), 
(2.21),  and  (2.22),  we  obtain 


dUih  d 


^  l\dx. 

axi  J^d 


[dr]  dr]  ^ 

[di ^ '‘'Wr 


J  x^=r) 


(  dd  ^ 

Ui  ^ - Uz 

V  dxi  J 


073  =  — d 


=0  [(2.21)] 


=0  [(2.20)] 


1  S')  , 


=o'](2.22)] 


(2.35) 


dUzh 

dt  dx 


fv  \  f  dr]  drj 

-J_^mU,dX3-^U,  +  - 


U3 


-I  Xs=r} 


(  dd  ^ 

U2  ^^1^ - '^S 

V  9xi  J 


J  xz^—d 


=0  [(2.21)] 


=0  [(2.20)] 


.  S')  , 

=  -h-r - h 

0X2 


(2.36) 


'^s2  —Tb2 

=01(2.22)) 

Rearranging  the  second  term  on  the  left  hand  sides  of  (2.35)  and  (2.36),  the  depth- 
integrated  horizontal  momentum  equations  are  given  by 


dU2h  _d 
dt  dxi  [ 


1U2  +  j  ^  (Wi  -  U\){U2  -  U2)  dxsj  =  “  '^*>2  (2-38) 


dt 
hUi 


with 


Ui  =  —  /  ui  dx3  =  depth-averaged  value  of  u\ 

h  J-d 

Ip 

C/2  =  —  «2  dxs  =  depth-averaged  value  of  1/2 

h  J-d 


(2.39) 

(2.40) 


In  the  following,  the  depth-integrated  continuity  equation  (2.33)  and  the  hor¬ 
izontal  momentum  equations  (2.37)  and  (2.38)  are  rewritten  using  the  conventional 
notations  of  x  =  xi,y  =  X2,u  =  «i,u  =  U2,U  =  Ui,  V  =  C/2,  and  r^y  = 

dh  d  . 


(2.41) 


(2.42) 

(2.43) 

m  =  /  (u  —  U^dz 

J 

(2.44) 

n  =  r  {u-  U){v  -  V)dz 

(2.45) 

-d 


The  dispersion  terms  m  and  n  defined  in  (2.44)  and  (2.45)  express  the  cross-shore 
and  alongshore  momentum  flux  corrections  due  to  the  vertical  variations  of  u  and 
V,  respectively.  The  dispersion  terms  result  from  the  vertical  integration  of  the 
horizontal  momentum  equations  {e.g.,  Rodi  1980).  It  is  noted  that  (2.41)  and  (2.42) 
with  m  =  0  are  the  same  with  those  used  previously  for  predicting  the  setup  and 
runup  of  normally  incident  waves  {e.g.,  Kobayashi  and  Wurjanto  1992). 

For  the  case  of  small  angles  of  incidence,  <  1,  the  dominant  cross-shore 
fluid  motion  computed  using  (2.41)  and  (2.42)  is  not  affected  by  the  secondary 
longshore  fluid  motion  varying  more  slowly  in  the  longshore  direction.  Furthermore, 
the  variations  in  the  j/-direction  appear  only  in  the  term  drifdy  in  (2.43)  and  along 
the  seaward  boundary  of  the  computation  domain.  The  alongshore  momentum 
equation  (2.43)  is  more  sensitive  to  the  gradual  alongshore  variability.  In  short,  the 
assumption  of  <C  1  reduces  computational  efforts  considerably  and  eliminates 
dilRculties  associated  with  lateral  boundary  conditions  for  general  two-dimensional 
computations. 

In  this  study,  the  numerical  model  is  developed  in  two  stages.  First,  the 
dispersion  terms  due  to  the  vertical  variations  of  the  horizontal  velocities  u  and  v 
are  neglected  by  setting  m  =  0  and  n  =  0.  This  2D  model  is  the  extension  of  the 
previous  one-dimensional  model  where  the  predicted  horizontal  velocities  are  the 
depth-averaged  velocities  only.  In  the  second  model,  the  dispersion  terms  m  and  n 
are  included.  The  vertical  profile  of  u  and  v  is  assumed  to  be  cubic  on  the  analogy 
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between  turbulent  bores  and  hydraulic  jumps.  For  brevity,  this  quasi-3D  model  is 
referred  to  as  the  ‘3D  model’. 

2.3  2D  Horizontal  Time-Dependent  Model 

As  mentioned  before,  (2.47)  and  (2.48)  with  m  =  0  are  the  same  as  the 
equations  used  previously  for  the  one  dimensional  model  of  Kobayashi  and  Wurjanto 
(1992).  They  stated  that  the  approximation  of  m  =  0  might  result  in  an  error  on  the 
order  of  10%.  Their  numerical  model  is  extended  herein  to  two  horizontal  dimensions 
using  the  approximation  n  =  0.  This  is  equivalent  to  neglecting  the  dispersion  terms 
arising  from  the  unknown  vertical  variations  of  the  horizontal  velocities.  The  depth- 
averaged  alongshore  velocity  V  is  computed  using  (2.43)  for  the  known  values  of  U 
and  h  computed  using  (2.41)  and  (2.42). 

For  the  2D  model,  the  physical  bottom  shear  stresses  and  r^y  are  assumed 
to  be  expressed  as 

rL  =  Ipfl'jwy  +  (V'Y  V  (2.46) 

=  \p}i^l{U'?  +  (V'r  V  (2.47) 

where  is  the  bottom  friction  factor  based  on  the  depth- averaged  velocities  {e.g., 
Kobayashi  et  al.  1987,  1989).  The  value  of  fl  specified  as  input  is  allowed  to  vary 
spatially  to  accommodate  the  spatial  variation  of  bottom  roughness  (Kobayashi  and 
Raichle  1994).  Normalization  of  the  bottom  shear  stresses  according  to  (2.5)  and 
(2.6)  with  (2.13)  yields 

n.  =  +  U  (2.48) 

n,  =  V  (2.49) 

For  <C  1,  the  normalized  bottom  shear  stresses  are  simplified  as 

T6.  =  /b|7/|C/  ;  ny  =  fi,\l7iv  (2.50) 
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with 


(2.51) 


2.4  Quasi-3D  Time-Dependent  Model 

In  this  quasi-3D  model,  the  dispersion  terms  m  and  n  are  included  to  account 
for  the  dispersion  effects  due  to  the  vertical  variation  of  the  horizontal  velocities 
u  and  V.  Two  new  equations  for  m  and  n  are  derived  from  the  corresponding 
three-dimensional  shallow-water  momentum  equations  (2.32)  using  the  algebraic 
procedure  based  on  the  method  of  moments. 

In  this  model,  the  bottom  boundary  layer  is  not  analyzed  explicitly  and  the 
normalized  bottom  stresses  for  <C  1  are  expressed  as 


nx  =  fb\uh\ub  ;  Tby  ^  fb\ub\vb  ;  fh  =  -^crf'b  (2.52) 

in  which  Ub  and  Vb  are  the  normalized  cross-shore  and  alongshore  velocities  immedi¬ 
ately  outside  the  bottom  boundary  layer,  respectively,  and  fl  is  the  bottom  friction 
factor  based  on  the  near-bottom  velocities  Ub  and  Vb  and  is  not  exactly  the  same  as 
the  bottom  friction  factor  fl  used  in  (2.46)  and  (2.47),  although  the  same  notation 
is  used  for  simplicity.  For  normally  incident  waves,  the  quadratic  friction  equation 
in  the  forms  of  (2.50)  and  (2.52)  was  shown  to  relate  the  bottom  shear  stress  and 
near-bottom  velocity  measured  inside  the  surf  zone  within  a  factor  of  two  (Cox  et 
al.  1996).  The  bottom  friction  factors  based  on  the  depth-averaged  velocity  and 
the  near-bottom  velocity  were  also  found  to  be  approximately  the  same  within  the 
accuracy  of  the  measurements  (Cox  et  al.  1996). 

To  derive  the  equation  for  m,  the  three-dimensional  momentum  equation 
(2.32)  with  i  =  I  multiplied  by  u  is  integrated  from  the  bottom  z  =  —d  to  the  free 
surface  z  =  r].  Applying  the  kinematic  boundary  conditions  (2.20)  and  (2.21)  at  the 
bottom  and  the  free  surface,  respectively,  yields 


d 


drj 


dr.. 
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The  integrated  equation  (2.53)  is  then  simplified  using  (2.41)  multiplied  by  to  ob¬ 
tain  the  depth-integrated  instantaneous  wave  energy  equation  derived  by  Kobayashi 
and  Wurjanto  (1992) 

=  p.54) 

The  specific  energy  E  defined  as  the  sum  of  kinetic  and  potential  energy  per  unit 
horizontal  area  is  given  by 

E  =  -  (?/^  +  +  m^  for  d  >  0 

E  =  -  (P  +  hU‘^  ^  m)  for  d<0  (2.55) 

in  which  the  potential  energy  is  taken  to  be  relative  to  the  potential  energy  in  the 
absence  of  wave  action  with  SWL  at  z  =  0.  The  energy  flux  Ep  per  unit  width  is 
expressed  as 

Ep  =  rjhU  +  ^{hU^  -h  dmU  -|-  m3)  (2.56) 

with  m3  =  kinetic  energy  flux  correction  due  to  the  third  moment  of  the  velocity 
deviation  {u~U)  over  the  depth  given  by 

m3=  [\u-Ufdz  (2.57) 

Jzb 

The  energy  dissipation  rate  D  per  unit  horizontal  area  in  (2.54)  is  given  by 

D=  f  T^^dz  (2.58) 

J-d  oz 

where  use  is  made  oi  =  0  ed  z  =  rj  and  no-slip  condition  «  =  0  at  z  =  —d.  With 
the  assumption  that  there  is  no  interaction  between  the  bottom  boundary  layer  and 
the  region  above  it,  the  energy  dissipation  rate  D  is  separated  into  Df  and  Dp 
(Kobayashi  and  Wurjanto  1992)  where  !>/  is  the  energy  dissipation  rate  inside  the 
bottom  boundary  layer  and  Dp  is  the  energy  dissipation  rate  outside  the  bottom 
boundary  layer  due  to  wave  breaking. 

i—d+Si  Qu  rv  du  , 

D  =  Df  +  Dp=  Ta;^dz+  '^x^dz  (2.59) 

J-d  oz  J-d-^St,  OZ 
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The  thickness  of  bottom  boundary  layer,  6b,  is  assumed  to  be  much  smaller  than 
the  water  depth  h=  (t)  +  d). 

It  has  been  mentioned  before  that  in  this  model  the  bottom  boundary  layer  is 
not  analyzed  explicitly.  The  energy  dissipation  rate  Dj  inside  the  bottom  boundary 
layer  may  be  approximated  by  (Jonsson  and  Carlsen  1976) 


j  ~~  ^bx'^h 


(2.60) 


By  use  of  (2.13),  the  energy  dissipation  rate  outside  the  bottom  boundary  layer  due 
to  wave  breaking  is  rewritten  as 


(2.61) 


where  the  lower  limit  z  =  —doi the  integration  should  be  interpreted  at  the  elevation 
immediately  outside  the  bottom  boundary  layer. 

Rearranging  the  instantaneous  wave  energy  equation  (2.54)  with  (2.55), 
(2.56),  (2.59)  and  (2.60),  and  by  use  of  the  depth-integrated  continuity  equation 
(2.41)  and  the  depth-integrated  momentum  equation  (2.42),  the  equation  for  the 
cross-shore  momentum  flux  correction  m  is  derived  as  (Kobayashi  and  Johnson 
1995) 

^  {ZmU  -t-  m3)  =  2  lu —  TbxUb  —  Db  j  (2.62) 

with 


Xib  'U/b  U 


(2.63) 


in  which  Ub  =  near  bottom  cross-shore  horizontal  velocity  correction  due  to  the  ver¬ 
tical  variation  of  the  cross-shore  horizontal  velocity  u  outside  the  bottom  boundary 
layer. 

To  derive  the  equation  for  n,  (2.32)  with  i  =  1  multiplied  by  v  and  (2.32)  with 
i  =  2  multiplied  by  u  are  added  and  integrated  from  z  =  —d  to  z  =  r}.  Applying  the 
kinematic  and  dynamic  boundary  conditions  at  the  free  surface  given  by  (2.20)  and 
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(2.22),  respectively,  and  the  kinematic  condition  (2.21)  and  the  no-slip  condition  at 
the  bottom  yields 

dr  ,  d  r  2..  t.rdv  ^rrdn  r(.dv^^du\, 


d  n 
dt  J-d 


r  uvdz+-f  r  uhdz = -hvp-  -  hup-  -  r 

J-d  dx  J-d  ox  oy  J-d 


which  can  be  rewritten  as 


^^{hUV  -f  n)  -h  ^{hU^V  +  2Un  +  Vm  +  na)  = 


,^,dri  ,  rrdri  r  f  dv  du\  , 

-  L  V’Vz  +  I 


(2.65) 


^3=  nu-Ufiv-V)^ 

J —d 


(2.66) 


where  m  and  n  are  defined  by  (2.44)  and  (2.45),  respectively. 

On  the  other  hand,  the  addition  of  (2.42)  multiplied  by  V,  (2.43)  multiplied 
by  U,  and  (2.41)  multiplied  by  (-UV)  yields 


Vnc,-UTky  (2.67) 


Subtracting  (2.67)  from  (2.65),  the  equation  for  the  alongshore  momentum  flux 

correction  n  is  obtained  as 

dn  d  ,  T/  ,  ^  jrd'm  ^,dn 

^  +  -^2Un  +  Vm  +  nz)-V—-U^ 


.r  rr  P  f  dv  du\  , 


(2.68) 


Similar  to  (2.59)  and  with  the  assumption  of  a  thin  boundary  layer,  the  last 
term  on  the  right  hand  side  of  (2.68)  may  be  approximated  as 


fv  (  dv  du\  r 


'I  /  dv  du\ 


(2.69) 


inside  thin  bottom  outside  bottom 

boundary  layer  boundary  layer 


By  use  of  (2.13),  (2.63)  and  (2.69),  (2.68)  can  be  rewritten  as 

dn  d  .  s.  rrdm  dU  _  .  _ 

~di'^  'dx  +  ’^3)  =  »/  -  n—  —  VlUx  -  UbThy  -  ^Dr, 


(2.70) 
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(2.71) 

(2.72) 


h 

Dn 


Vb 


-I 


V 
du  dv 


-d  az  oz 


in  which  Vh  =  alongshore  horizontal  velocity  correction  due  to  the  vertical  variation 
of  the  alongshore  horizontal  velocity  v  immediately  outside  the  bottom  boundary 
layer. 

Since  the  thickness  of  the  bottom  boundary  layer  is  assumed  to  be  much 
smaller  than  the  water  depth  /i  =  (?/  +  d),  the  lower  limit  z  =  — d  of  the  integration 
in  (2.61)  and  (2.72)  should  be  interpreted  at  the  elevation  immediately  outside  the 
bottom  boundary  layer.  The  contributions  of  the  boundary  layer  flow,  which  must 
satisfy  u  =  0  and  u  =  0  at  the  bed,  to  the  second  moments  m  and  n  given  by 
(2.44)  and  (2.45)  and  the  third  moments  m3  and  nz  defined  by  (2.57)  and  (2.66), 
respectively,  are  assumed  to  be  negligible. 

To  obtain  h,U,V,m  and  n  using  (2.41)  -  (2.43),  (2.62)  and  (2.70),  the  vari¬ 
ables  U6,  Uf,,  m3,  ^3,  Db  and  need  to  be  expressed  in  terms  of  these  five  unknown 
variables.  As  a  first  attempt  to  deal  with  this  closure  problem,  the  instantaneous 
horizontal  velocities  u  and  v  outside  the  bottom  boundary  layer  are  assumed  to  be 
expressed  as 

u{t,x,y,z)  =  U{t,x,y)  +  Ui{t,x,y)F{C)  (2.73) 


and 

v{t,x,y,z)  =  V{t,x,y)  +  Vb{t,x,y)F{C)  (2.74) 


with 


^  _  z  +  d{x,y) 
h{t,x,y) 


(2.75) 


in  which  F  is  assumed  to  be  a  function  of  (  only,  with  C  =  0  at  z  =  —d  and  C  =  1 
at  z  =  rj.  The  definitions  of  Ui  and  Vb  in  (2.63)  and  (2.71)  require  F  =  1  at  the 


bottom  C  =  0-  Furthermore,  the  dimensional  turbulent  eddy  viscosity  v[  outside  the 
bottom  boundary  layer  is  assumed  to  be  given  by 


/^2 


i/,'  =  {C,h') 


du' 


dz' 


(2.76) 


in  which  Cn  —  mixing  length  parameter.  The  turbulence  measurements  inside  the 
surf  zone  by  Cox  et  al.  (1994)  indicate  that  the  mixing  length  parameter  Ct  is 
on  the  order  of  0.1.  Accordingly,  the  eddy  viscosity  vt  normalized  in  (2.6)  may  be 
expressed  as 
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Vt  =  crCih 


du 


dz 


(2.77) 


Substitution  of  (2.73),  (2.74),  and  (2.77)  into  (2.44),  (2.45),  (2.57),  (2.61),  (2.66) 
and  (2.72)  yields 


m  =  C2hul 

,  n  —  j 

(72=  f^F^dC 

Jo 

(2.78) 

mz  =  Czhul 

;  nz  =  CzhubVb  ; 

CO 

II 

(2.79) 

Db  =  CB<^Cj\ub  1^  ; 

Dn  =  CbctC^  I  Ub\ubh 

II 

dF 

dC 

3 

dC 

(2.80) 

in  which  the  constants 

,  Cz  and  Cb  can  be  found  for  the  specified  shape 

function 

F.  To  find  Ub  using  (2.78)  for  the  computed  h>  0  and  m  >  0,  it  is  assumed  that 


Ub  •“ 

/  m  \  2 

[CzhJ 

for 

U>0 

Ub  = 

/  m  \2 

KCzh) 

for 

U  <0 

(2.81) 

to  ensure  that  |«i,|  <  \U\  where  Ub  —  {U  +  U6)  is  the  near-bottom  cross-shore  velocity 
used  in  (2.52).  After  uj,  is  obtained,  then  Vb  =  nl{C2hub)^  Vb  =  {V  +  Vb)i  and  (2.79) 
and  (2.80)  yield  m3,n2,DB  and  £)„. 

Finally,  the  shape  function  F  needs  to  be  specified.  Svendsen  and  Madsen 
(1984)  assumed  a  cubic  velocity  profile  for  their  analysis  of  a  single  turbulent  bore 


25 


on  a  beach.  For  regular  and  irregular  breaking  waves  on  beaches,  the  following  cubic 
profile  is  tentatively  assumed: 

F  =  1  -  (3  +  0.75a)  +  aC®  for  0  <  C  <  1  (2.82) 


in  which  a  =  cubic  velocity  profile  parameter.  Equations  (2.73)  and  (2.74)  require 
that  F  =  1  at  (  =  0  and  the  depth-averaged  value  of  F  equals  zero.  Substitution 
of  (2.73)  -  (2.75)  and  (2.77)  into  (2.13)  yields  the  shear  stresses  Tx  and  Ty.  Equa¬ 
tion  (2.82)  predicts  zero  shear  stresses  immediately  outside  the  bottom  boundary 
layer  in  contradiction  with  the  turbulence  measurements  inside  the  surf  zone  by  Cox 
et  al.  (1994).  Moreover,  the  shear  stresses  at  the  surface  are  zero  only  if  a  =  4. 
Consequently,  (2.82)  with  the  single  parameter  a  may  not  predict  the  shear  stresses 
accurately  in  the  vicinity  of  the  free  surface  and  bottom.  Comparison  of  (2.82)  and 
the  cubic  profile  assumed  by  Svendsen  and  Madsen  (1984)  suggests  that  a  is  about  3. 

Substitution  of  (2.82)  into  (2.78)  -  (2.80)  yields  the  explicit  expressions  of 
C2,  Cz  and  Cb  in  terms  of  a  (Kobaycushi  and  Johnson  1995). 
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a 


C2  =  l  +  '7r  +  7:  + 


6^  ah 


C3 


Cb 


5  3  7 

3a  36^  ,  3a^  -f  6^  3a6^  aH 

l  +  5+-H--  +  a6-K-^  +  — +  -+^ 

/  ,,  36a62  ^  2,  27a^\ 

-{2b^  +  -—  +  9aH+  —  ] 


where 


b  =  -(3 -h  0.75a) 


(2.83) 

(2.84) 

(2.85) 

(2.86) 


For  the  range  of  a  =  3  -  4,  the  assumed  cubic  profile  F  is  not  very  sensitive  to  the 
parameter  a.  Furthermore,  C2  =  0.49  -  0.55,  Cz  =  -0.07  -  0.00,  and  Cb  =  12.3  -  15.2 
for  a  =  3  -  4.  The  computed  results  using  Ci  =  0.1  -  0.2  in  (2.80)  and  a  =  3  -  4  in 
(2.83)-(2.85)  are  found  to  be  very  similar  (Johnson  et  al.  1996).  The  typical  values 
of  Ce  =  0.1  and  a  =  3  are  hence  employed  for  the  3D  computed  results  presented 
in  Chapters  5  and  6. 
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2.5  Time- Averaged  Continuity  and  Momentum  Equations 

The  time-dependent  models  based  on  the  depth-integrated  continuity  equa¬ 
tion  (2.41),  the  depth-integrated  cross-shore  momentum  equation  (2.42),  and  the 
depth-integrated  alongshore  momentum  equation  (2.43)  allow  us  to  compute  the 
oscillatory  and  mean  components  of  the  water  depth  /i,  depth-averaged  cross-shore 
velocity  t/,  and  depth-averaged  alongshore  velocity  V.  The  time-averaged  quantity 
denoted  by  an  overbar  is  defined  as 

M  = - - -  M{t)  dt  (2.87) 

^end  ^begin  *begiii 

where 

M  =  computed  time- varying  quantity  at  given  location 
^begin  =  normalized  time  when  the  time  averaging  begins 
^end  =  normalized  time  when  the  time  averaging  ends 

The  time-averaged  continuity  equation  derived  from  (2.41)  is  given  by 

hU  =  0  (2.88) 


which  satisfies  the  condition  of  no  flux  into  the  impermeable  beach.  The  time- 
averaged  cross-shore  momentum  equation  derived  from  (2.42)  with  the  assumption 
that  m  =  0  and  /j  =  0  is  similar  to  the  conventional  equation  used  to  predict  wave 
setup  (Kobayashi  et  ah  1989).  The  time-averaged  alongshore  momentum  equation 
corresponding  to  (2.43)  can  be  written  as 


d  _  dn  _  j-dri  1  d- - 

dx  ^dy  2 


(2.89) 


with 


=  hUV 


(2.90) 


If  the  mean  and  variance  of  the  free  surface  elevation  y  above  SWL  do  not  vary  in 
the  alongshore  direction,  (2.89)  is  similar  to  the  conventional  alongshore  momentum 


equation  in  which  ~  0  is  assumed  and  the  dispersion  term  n  is  expressed  in  the 
form  of  lateral  mixing  {e.g.,  Longuet-Higgins  1970).  In  this  study,  n  =  0  for  the  2D 
model  and  n  is  explicitly  computed  using  (2.70)  for  the  3D  model.  The  term  Sxy 
given  by  (2.90)  may  be  regarded  as  the  alongshore  radiation  stress  (e.g.,  Thornton 
and  Guza  1986)  where  Sxy  —  h(U  —  U){V  —  V)  for  7/  ~  0  and  hU  =  0. 


2.6  Time- Averaged  Energy  Equation 

The  time-averaged  cross-shore  energy  equation  corresponding  to  (2.54)  and 
(2.59)  is  given  by 

^  (2-91) 

ax 


with 


AE  = 


E{t  —  tend)  E{t  —  tbegin) 


^end  ^b' 


egin 


(2.92) 


For  regular  waves,  AE  defined  by  (2.92)  is  zero  because  of  periodicity,  whereas  for 
irregular  waves,  tbegin  and  tend  are  chosen  such  that  AE  is  negligible  to  obtain  the 
time-averaged  energy  equation  for  a  stationary  sea  state.  As  a  result,  AE  can  be 
deleted  from  (2.91)  but  is  included  for  completeness. 

In  the  previous  one-dimensional  model  and  the  present  2D  model,  Db  is 
calculated  using  (2.91)  because  these  models  do  not  include  any  physical  dissipa¬ 
tion  mechanism  associated  with  wave  breaking.  For  the  3D  model,  the  values  of 
E,EF,Df  and  Db  are  computed  using  (2.55),  (2.56),  (2.60)  and  (2.80),  respec¬ 
tively.  The  time-averaged  energy  equation  (2.91)  will  be  used  to  check  the  degree 
of  numerical  dissipation  in  Chapters  5  and  6. 

Equations  (2.61)  and  (2.62)  give  Db  =  0  if  u  is  independent  of  z  so  that 
u  =  =  0,m  =  0,  and  m3  =  0.  This  proves  that  the  energy  dissipation 

due  to  wave  breaking  in  the  previous  one-dimensional  model  and  the  present  2D 
model  based  on  the  assumptions  of  Ui,  =  0,  m  =  0,  and  m3  =  0  is  solely  numerical 
(Kobayashi  and  Wurjanto  1992).  The  normalized  energy  dissipation  rate  Db  due 
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to  the  vertical  variations  of  and  u  outside  the  boundary  layer  as  given  by  (2.61) 
is  the  same  as  the  dissipation  rate  due  to  breaking  of  normally  incident  waves  used 
by  Svendsen  and  Madsen  (1984). 

In  addition  to  the  nonlinear  shallow  water  equations  used  in  this  study,  the 
motion  of  waves  in  the  nearshore  region  is  also  modeled  using  the  Boussinesq  equa¬ 
tions.  To  include  energy  dissipation  due  to  wave  breaking  in  the  Boussinesq  equa¬ 
tions,  Zelt  (1991)  added  a  term  corresponding  to  the  term  m  in  (2.42)  and  expressed 
this  term  in  the  form  of  horizontal  momentum  diffusion  with  an  artificial  viscosity 
that  was  calibrated  empirically  for  breaking  solitary  waves.  Schaffer  et  al.  (1992)  ex¬ 
pressed  the  additional  momentum  fluxes  similar  to  the  terms  m  and  n  in  (2.42)  and 
(2.43)  using  a  simple  surface  roller  model  in  which  an  empirical  geometric  method 
was  used  to  estimate  the  shape  and  location  of  the  surface  roller.  However,  Zelt 
(1991)  and  Schaffer  et  al.992)  did  not  check  whether  their  computed  results  satis¬ 
fied  the  wave  energy  equation.  Consequently,  the  degree  of  numerical  dissipation  in 
their  computed  results  was  not  certain. 
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Chapter  3 


NUMERICAL  MODELS 


The  governing  equations  for  the  2D  and  3D  models  are  solved  numerically 
using  two  different  numerical  methods.  The  unknown  variables  for  the  cross-shore 
motion  in  the  2D  model  are  computed  using  the  dissipative  Lax-Wendroff  method 
in  the  same  way  as  the  one-dimensional  model  explained  by  Kobayashi  et  al.  (1987, 
1989)  since  the  governing  equations  for  the  cross-shore  motion  for  the  2D  model  are 
the  same  as  those  of  the  one-dimensional  model.  The  alongshore  motion  in  the  2D 
model  as  well  as  the  cross-shore  and  alongshore  motions  in  the  3D  model  are  solved 
using  the  MacCormack  method  (MacCormack  1969).  This  method  is  a  simplified 
variation  of  the  two-step  Lax-Wendroff  method  {e.g.,  Anderson  et  al.  1984)  and 
has  been  applied  successfully  for  the  computation  of  unsteady  open  channel  flows 
with  hydraulic  jumps  {^e.g.,  Fennemaand  Chaudhry  1986;  Gharangik  and  Chaudhry 
1991). 

The  two  models  use  the  same  coordinate  system,  boundary  and  initial  condi¬ 
tions  as  well  as  the  specified  incident  wave  train  at  the  seaward  boundary,  although 
the  3D  model  requires  the  additional  conditions  for  the  momentum  flux  correction 
terms  m  and  n.  Figure  3.1  shows  the  finite  difference  grid  of  constant  grid  sizes  Aa: 
and  Aj/  used  to  solve  the  governing  equations  (2.41),  (2.42)  and  (2.43)  numerically. 
The  cross-shore  coordinate  x  is  taken  to  be  positive  landward.  The  alongshore  coor¬ 
dinate  y  is  positive  in  the  downwave  direction  with  y  =  0  at  the  upwave  boundary. 
The  cross-shore  lines  in  Figure  3.1  are  located  at  y  =  (i  —  l)Ay  with  i  =  1,  2,  . . . ,  7 
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where  I  is  the  number  of  the  cross-shore  lines.  The  seaward  boundary  of  the  compu¬ 
tation  domain  is  located  at  a;  =  0  along  the  y-axis.  Along  each  cross-shore  line,  the 
nodes  are  located  at  a:  =  {j  —  l)Aa:  with  j  =  1,  2,  . . . ,  J  where  J  is  the  number  of 
the  cross-shore  nodes.  The  computational  shoreline  is  defined  as  the  location  where 
the  normalized  instantaneous  water  depth  h  equals  a  small  value  S  as  in  the  previ¬ 
ous  one-dimensional  model  .  The  initial  time  t  =  0  for  the  computation  marching 
forward  in  time  is  taken  to  be  the  time  when  the  incident  wave  train  arrives  at  x  =  0 
and  there  is  no  wave  action  in  the  region  x  >  0  and  y  >  0. 


\  incident 
wave 


y 


moving 

shoreline 


Figure  3.1:  Finite  difference  grid  for  the  numerical  models. 


The  incident  wave  free  surface  elevation,  rii{t,y),  at  x  =  0  and  y  =  (i  —  l)Ay 
with  i  =  1,  2,  . . . ,  /  needs  to  be  specified  as  input  to  the  numerical  model.  The 
bottom  profile  along  each  of  the  /  cross-shore  lines  needs  to  be  specified  as  well 
unless  the  bottom  profile  is  assumed  to  be  uniform  in  the  alongshore  direction.  The 
value  of  d  =  {h  —  ri)  is  the  normalized  depth  below  SWL  and  is  known  for  the 
specified  bottom  profile. 

In  Chapters  4,  5  and  6,  the  2D  and  3D  numerical  models  will  be  compared 
with  available  data  whose  beach  profile  and  incident  wave  conditions  may  be  as¬ 
sumed  to  be  uniform  in  the  alongshore  direction.  For  this  case,  it  is  sufficient  to 
compute  the  unknown  variables  in  the  cross-shore  direction  along  three  lines  at 
y  =  0,  Ay  and  2Ay  with  7  =  3  in  Figure  3.1,  and  then  compute  the  unknown  vari¬ 
ables  involved  in  the  alongshore  wave  motion  along  the  center  line  at  y  =  Ay.  The 
computed  y  along  the  three  lines  are  used  to  ensure  the  alongshore  uniformity  of 
the  mean  and  variance  of  y  appearing  in  the  time-averaged  alongshore  momentum 
equation  (2.89). 

It  is  also  noted  that  even  if  the  beach  profile  and  incident  wave  conditions 
vary  gradually  alongshore,  the  developed  numerical  models  should  be  applicable  for 
the  computation  of  the  gradual  longshore  variations  of  the  wave  motion  by  choosing 
a  larger  value  of  7,  provided  that  lateral  boundary  conditions  do  not  affect  the  wave 
motion  in  the  computation  domain.  However,  computations  for  this  more  general 
case  have  not  been  made  yet  and  beyond  the  scope  of  this  study. 

3.1  2D  Numerical  Model 

3.1.1  Summary  of  Equations 

The  equations  used  for  the  2D  model  are  summarized  in  the  following.  The 
depth-integrated  continuity  equation  (2.41)  and  the  horizontal  momentum  equations 
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(2.42)  with  m  =  0  and  (2.43)  with  n  =  0  are  rewritten  as 

0 


'^bx 


dri 

^dy 


(3.1) 

(3.2) 

(3.3) 


with 

n.  =  /i|£i|£i  1  n,  =  m\V  1  A  =  i<7/;  (3.4) 

The  time-averaged  alongshore  momentum  equation  (2.89)  with  n  =  0  is  given 


by 


with 


-h^ 

dx^^~  ny  riQy 


1  d-, - z:;! 

2 


(3.5) 


S^y  =  hUV 


(3.6) 


where  overbar  denotes  time-averaged  quantities  defined  in  (2.87). 

The  time-averaged  cross-shore  energy  equation  (2.91)  is  given  by 

AE  +  ^(1F)  = -:d7-:D^  (3.7) 

with 


E  = 

i  (t/^  +  hU^)  for  d  >  0 

E  = 

i  (t)^  -<f  +  hU^)  for  d  <  0 

(3.8) 

AE  = 

E{t  =  tend)  -  E{t  =  tbegin) 

(3.9) 

^end  ^begin 

Ejr  = 

rjhU  +  ^{hU^) 

(3.10) 

Df  = 

UxU 

(3.11) 

Db  is  calculated  using  (3.7). 


3.1.2  Numerical  Method  for  Cross-Shore  Wave  Motion 

For  the  known  values  of  h,r],U  and  V  at  the  time  level  t  and  at  all  the 
nodes  used  in  the  computation,  the  values  of  these  variables  at  the  next  time  level 
t*  =  (t  +  At),  which  is  denoted  by  the  superscript  asterisk,  are  computed  in  se¬ 
quence.  The  time  step  size  At  is  allowed  to  vary  in  a  manner  similar  to  the  existing 
one-dimensional  model  where  At  is  reduced  in  a  semi-automated  way  whenever 
numerical  difficulties  occur  at  the  moving  shoreline. 

Since  the  governing  equations  of  the  2D  model  for  the  cross-shore  wave  motion 
with  small  angles  of  wave  incidence  are  the  same  as  those  of  the  one-dimensional 
model  along  each  of  the  cross-shore  lines,  the  values  of  h'j  and  Uj  at  the  node  j  with 
j  =  1,2,  s*  are  computed  from  (3.1)  and  (3.2)  using  the  dissipative  Lax-Wendroff 
method  in  the  same  way  as  the  one- dimensional  model  explained  by  Kobayashi  et  al. 
(1987,  1989).  The  integer  s*  indicates  the  wet  node  next  to  the  moving  shoreline  at 
the  next  time  level  t*.  The  incident  wave  train  rji  specified  at  x  =  0  for  each  of  the 
cross-shore  lines  is  the  input  to  this  cross-shore  fluid  motion  computation. 

The  Lax-Wendroff  finite  difference  method  is  discussed  briefly  in  the  following 
for  the  subsequent  comparison  between  this  method  and  the  MacCormack  method 
used  to  solve  the  alongshore  momentum  equation  (3.3).  The  continuity  and  cross- 
shore  momentum  equations  (3.1)  and  (3.2)  and  (3.4)  can  be  rewritten  in  terms  of  h 
and  q  =  hU  as  follows: 


dq 

dt  dx 


dh  dq 
dt  dx 


0 

-eh-h\u\u 


(3.12) 

(3.13) 


where  9  =  d{—d)/dx  is  the  normalized  local  slope. 

Equations  (3.12)  and  (3.13)  can  be  combined  and  expressed  in  the  following 
vector  form: 


dV  dF  ^  ^ 


(3.14) 
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where 


U 

F 

G 


9 


Gi 

G2 


eh  +  fi\u\u 

0 


(3.15) 

(3.16) 

(3.17) 


Using  a  Taylor  series  expansion,  U*  at  the  next  time  level  t*  =  {t  +  At)  is  expressed 
as 

V*{x,t  +  At)  =  U(x,t)  +  ^  ^ 


The  main  feature  of  the  Lax-Wendroff  scheme  of  second-order  accuracy  in 
At  and  Ax  is  to  express  the  temporal  derivatives  in  (3.18)  in  terms  of  the  spatial 
derivatives.  From  (3.14),  the  temporal  derivatives  of  U  can  be  written  as 

^  =  (3.19) 

dt  dx  ^  ^ 


and 


dt^  dt  \  dx 


in  which  dFfdt  can  be  expressed  as 


dx\dt)  dt 


(3.20) 


^ 

dt  dt 


(3.21) 


where 


(3.22) 


By  use  of  (3.12)  and  (3.13),  dGildt  can  be  shown  to  be  expressed  as 


dGi 

dt 


1£1 

h 


(3.23) 
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Substitution  of  (3.19)  and  (3.20)  together  with  (3.21)  into  (3.18)  leads  to 

+  0{Atf  (3.24) 


• 

’aF 

At  2 

■  d 

dG 

U*  =  U  -  At 

+  ^ 

2 

dx 

dt 

Using  the  finite  difference  approximations  of  second-order  accuracy  in  Ax,  (3.24) 
can  be  shown  to  be  expressed  in  the  following  explicit  form: 


U*  =  U,-  -  A 


-  Fj_i)  -1-  AxGj 


+  -  g,.,  -  AxSi)  +  D,-  (3.25) 


with 


A  = 


At 

Ax 


1  r  /\x 

gj  =  +  A,)  -  F,- +  — (G,+i  +  G,-) 

Axej  -  -  qj-i) 


(3.26) 

(3.27) 

(3.28) 


ej  =  2(fi)j 


hj 


(uf  -  h,y 


■Uj 


9j+i  ~  gj-i 
2Ax 


-  -  imujiuj 


(3.29) 


The  superscript  asterisk  indicates  the  quantities  at  the  next  time  level  t*  =  {t  +  At) 
and  the  subscript  j  denotes  the  nodal  location  at  Xj  =  {j  —  l)Ax. 

The  vector  in  (3.25)  represents  the  additional  term  for  damping  high 
frequency  parasitic  oscillations,  which  tend  to  appear  at  the  rear  of  a  breaking 
wave,  and  is  given  by  (Kobayashi  et  al.  1987) 


with 


D,  =  5  [Qi(Ui+i  -  U>)  -  -  U,_,)] 


Qi  =  Pjl  + 


(3.30) 


(3.31) 


where  I  =  unit  matrix.  The  coefficients  pj  and  Vj  are  given  by 
1 


Pi  = 


2(cj+i  +  Cj) 

£i|^j-n  ^214’ j+i  -  i’j 

{Cj+1  +  Cj) 


[ei\4>j+i  -  (j>j\{ipj+i  +  ‘4’j)  -  e2l'0j+i  -  V’j|(<^j+i  +  (3-32) 

(3.33) 
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with 


c=y/h  ;  <j)=U  +  c  ;  tj;  =  U  —  c 

where  ei  and  €2  are  the  positive  numerical  damping  coefficients  on  the  order  of  unity. 


3.1.3  Numerical  Method  for  Alongshore  Velocity 

Along  each  of  the  cross-shore  lines  excluding  the  cross-shore  lines  ai  i  =  1 
and  J,  the  value  of  V*  at  the  node  j  with  j  =  1,  2,  . . . ,  s*  is  computed  using  (3.3), 
which  is  expressed  in  terms  of  qi  =  hV  where  =  hjVj  and  q}.  =  h*Vj 


This  equation  for  qi,  is  solved  numerically  using  the  MacCormack  method  (MacCor- 
mack  1969).  The  use  of  the  MacCormack  method  eliminates  the  algebraic  manipu¬ 
lations  involved  in  the  Lax-Wendroff  method  to  express  the  temporal  derivatives  in 
terms  of  the  spatial  derivatives  performed  in  (3.19),  (3.20)  and  (3.21). 

The  MacCormack  method  consists  of  two  steps,  z.e.,  predictor  and  correc¬ 
tor  steps,  with  forward  and  backward  spatial  differences  which  are  interchangeable 
(Anderson  et  al.  1984).  When  applied  to  (3.34),  the  predictor  step  is  expressed  as 


=  %  -  £  -  u’%)  -  Ai 


'•i 


,  dy 


(3.35) 


for  J  =  1,  2,  . . . ,  s*  where  C/J  and  h*  have  already  been  computed  and  =  0 

landward  of  the  moving  shoreline  node  s* .  The  term  qt-  given  by  (3.35)  is  a  tem¬ 
porary  predicted  value  at  the  next  time  level  t*.  The  normalized  bottom  friction 
factor  fb  given  by  (3.4)  is  allowed  to  vary  spatially  and  {fb)j  is  used  in  (3.35).  The 
corrector  step  of  the  MacCormack  method  is  given  by 


% 


At 
~  ~Ax 


'dri* 

.  9y 


+ (Mi \m‘.  ('>;) 


-1 


(3.36) 
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for  j  =  2,  3,  . . . ,  5*.  Finally,  the  value  of  at  the  node  j  and  at  the  next  time 
level  t*  is  computed  from 


9^. 


j  =  2,3,...,  s* 


with 


2Ax 


[ic'rt.  - c/;i  {%„  - «,)  -  \U-  - u-_,\  („,  - 


(3.37) 


(3.38) 


where  Dj  is  the  numerical  damping  term  based  on  the  procedure  described  in  Richt- 
myer  and  Morton  (1967)  and  is  the  numerical  damping  coefficient.  The  term  Dj 
is  added  to  (3.37)  except  for  j  =  s*  and  =  0  with  j  =  (s*  +  1),  . . . ,  J  landward 
of  the  moving  shoreline. 

The  term  {dr)*/dy)j  for  the  cross-shore  line  i  in  (3.35)  and  (3.36)  is  approxi¬ 
mated  by  a  central  finite  difference  based  on  the  computed  values  of  rjj  at  the  two 
adjacent  cross-shore  lines  (f  —  1)  and  {i  -f  1)  if  all  the  nodes  j  at  the  three  lines 
{i  —  l),f  and  (i  +  1)  are  wet  and  seaward  of  the  shoreline  whose  location  varies  for 
the  three  cross-shore  lines.  Otherwise,  this  term  is  set  to  be  zero  because  the  along¬ 
shore  fluid  motion  turns  out  to  be  sensitive  to  spurious  shoreline  oscillations  of  a 
short  duration  unlike  the  cross-shore  fluid  motion  as  will  be  explained  in  Section  4.3 
in  relation  to  the  computed  results  of  V  and  U. 

The  numerical  damping  term  Dj  given  by  (3.38)  tends  to  reduce  high-frequen¬ 
cy  numerical  oscillations  at  the  rear  of  a  steep  wave  front  without  modifying  the 
slowly  varying  part  of  the  alongshore  fluid  motion.  The  computed  alongshore  fluid 
motions  with  the  numerical  damping  coefficient  =  0  and  1  are  found  to  be  the 
same  except  for  the  slight  reduction  of  the  high-frequency  numerical  oscillations  as 
will  be  shown  in  Section  4.3.  The  computed  results  presented  in  Sections  4.2  and 
4.3  are  based  on  =  1  unless  stated  otherwise. 

It  is  noted  that  a  forward  spatial  difference  is  used  for  the  second  term  on  the 
right  hand  side  of  the  predictor  equation  (3.35),  while  a  backward  spatial  difference 
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is  used  for  this  term  in  the  corrector  equation  (3.36).  This  is  because  the  spatial 
diiference  in  the  predictor  equation  is  recommended  to  be  in  the  direction  of  prop¬ 
agation  of  wave  fronts  (Anderson  et  al.  1984).  Since  wave  breaking  can  also  occur 
during  wave  downrush  on  the  slope,  the  computation  is  also  made  by  reversing  the 
spatial  differencing  in  the  predictor  and  corrector  equations.  The  computed  results 
are  found  to  be  practically  the  same  as  will  be  shown  in  Section  4.3.  As  a  result, 
the  MacCormack  method  based  on  (3.35)  and  (3.36)  should  be  applicable  to  wave 
fronts  propagating  seaward  as  well  as  landward. 


3.1.4  Numerical  Stability  Criterion  for  MacCormack  Method 

The  linear  stability  criterion  for  the  MacCormack  method  applied  to  (3.34) 
without  the  terms  on  the  right  hand  side  of  (3.34)  is  the  Courant  condition  {e.g., 
Ryrie  1983;  Anderson  et  al.  1984)  as  shown  in  the  following.  For  the  linear  stability 
criterion,  U  is  assumed  constant.  The  predictor,  corrector  and  final  steps  of  the 
MacCormack  method  applied  to  this  simplified  alongshore  momentum  equation  may 
be  written  as 

UAt  I  X 

qti  =  %  -  -  qej) 

UAt.  .  N 

qt-j  qtj  q^j—i ) 

ql  =  \{qi,  +  %) 

Substituting  ije,.  into  the  final  step,  the  value  of  is  given  by 
qij  =  2  ~ 


with 


cx  = 


UAt 

Ax 
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For  the  assumed  form  of  with  Xj  =  {j  -  l)Aa:,  the  amplification  factor 

G  for  the  time  step  from  the  present  time  t  to  the  next  time  level  t*  =  {t  +  At)  is 
written  as 

G  =  =  1  —  a^{l  —  cos  kAx)  —  ia  sin  kAx 

The  stability  condition  of  |G|  <  1  is  satisfied  if  a  <  1  (Anderson  et  al.  1984), 
which  leads  to 

\Ura\^t 


Ax 


<  1 


(3.39) 


where  \Um\  is  the  maximum  absolute  value  of  U  in  the  computation  domain. 

On  the  other  hand,  the  numerical  stability  criterion  for  the  dissipative  Lax- 
Wendroff  method  used  to  solve  (3.12)  and  (3.13)  is  given  by  (Packwood  1980) 

(3.40) 


At 


< 


Ax  T  c„ 


with 


Cyn  =  maximum  value  of  y/h  >  0 

e  =  maximum  value  of  e\  and  £2  introduced  in  (3.32)  and  (3.33). 
Using  the  following  inequalities 

2  .1 

<  1  for  e  >  0 


and 


(3.40)  leads  to 


\Um\  +  Cm  >  \Urt 


At  1 


(3.41) 


As  a  result,  the  Courant  condition  for  the  MacCormack  method  used  to  solve 
(3.37)  is  satisfied  if  the  numerical  stability  criterion  for  the  dissipative  Lax-Wendroff 
method  used  to  solve  (3.12)  and  (3.13)  is  satisfied. 

For  the  computed  results  presented  in  Chapter  4,  the  value  of  At  smaller 
than  that  required  by  the  numerical  stability  criterion  is  used  to  minimize  spurious 
shoreline  oscillations  as  will  be  discussed  in  more  detail  in  Section  4.3. 
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3.1.5  Seaward  Boundary  Algorithm 

The  value  of  with  j  =  1  at  the  seaward  boundary  a;  =  0  can  not  be 
computed  using  (3.37).  To  devise  an  appropriate  seaward  boundary  algorithm,  the 
depth-integrated  alongshore  momentum  equation  (3.34)  is  expressed  in  the  following 
characteristic  form  as  depicted  in  Figure  3.2: 


dv  dv  _  dri  Mu\v 

dt  dx  dy  h 


(3.42) 


where  use  is  made  of  the  continuity  equation  (3.12).  The  variation  of  the  character¬ 
istic  variable  V  along  dxfdt  =  U\s  given  by  (3.42). 

To  compute  Vj*  and  =  hlV^  using  (3.42),  the  linear  approximation  of 
dxjdt  =  U  at  X  =  0  and  at  the  time  level  T  is  expressed  as  SxjAt  =  as  depicted 
in  Figure  3.2  where  the  spatial  increment  6x  satisfies  |dx|  <  Ax  due  to  the  Courant 
condition,  At  <  Ax/\Um\-  If  >0,  dx  >  0  and  the  characteristic  path  is  directed 
landward.  If  <0,  dx  <  0  and  the  characteristic  path  is  directed  seaward. 
The  finite  difference  approximation  of  (3.42)  along  this  characteristic  path  may  be 
expressed  as 

T/*  __  t7  I'f,  ^  I  TV*  IT/* 

(3.43) 


V{-V 


'dv* 


At 


dy 


K 


where  V  is  the  value  of  V  at  the  start  of  this  characteristic  path  at  the  time  level  t 
with  dx  =  At  Ui- 

If  >  0,  the  characteristic  variable  propagates  into  the  computation  domain 
and  V  should  be  specified  as  input;  however,  this  option  is  not  feasible  for  practical 
applications.  Alternatively,  V  may  be  estimated  from  Vi  and  V2  by  interpolation 
for  <  0  and  by  extrapolation  for  >  0. 


V  =  Vi-  {V2  -  Vr) 


Substitution  of  (3.44)  into  (3.43)  yields 

n  -1  r 


Vi*  = 


1  + At 


iJMS' 


k; 


(3.44) 


(3.45) 
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Figure  3.2:  Sketch  of  seaward  boundary  algorithm  for  alongshore  velocity  V. 

where  {dTj*/dy)i  is  approximated  by  a  central  finite  difference  based  on  the  computed 
values  T/i  at  the  two  adjacent  cross-shore  lines. 

3.2  3D  Numerical  Model 
3.2.1  Summary  of  Equations 

The  equations  used  for  the  quasi-3D  model  are  summarized  in  the  follow¬ 
ing.  The  depth-integrated  continuity  equation  (2.41)  and  the  horizontal  momentum 
equations  (2.42)  and(2.43)  are  given  by 


• 

dh  d  ,,  rrN 

=  0 

-  ^  ^hx 

OX 

-  ^Qy 

(3.46) 

(3.47) 

(3.48) 

• 

with 

m  =  J  {u  —  Uydz  ;  =  J 

(u  -  U){v  -  V)dz 

—d 

(3.49) 

'^hx  —  7  '^hy  j  fb  2,^^^ 

(3.50) 

ui,  =  U  +  ub  ;  Vi 

=  V  +  Vb 

(3.51) 
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The  two  additional  equations  (2.62)  and  (2.70)  for  m  and  n  are  given  by 

(3.52) 


—  +  —  (3mt/  +  m3)  =  2(7/  —  -  UbUx  -  Db 


with 


dn  d.  .  .  ,,dm  dU  _  . 

^  [nil  +  ml/  +n3)  =  V  —  -  n—  -  VbUx  -  UbUy  -  2L>„ 


=  -  1  7TT  I '  for  7/  >  0  ;  ut  =  ( 7^ ) "  for  7/  <  0  ; 


(3.53) 


7/2/1 


n 


Vb  = 


C2hub 


;  7/2=  C  FUC 

Jo 


mz  =  Czhul  ;  nz  =  Czhulvb  ;  Cz  =  f  F^  d( 

J  0 


(3.54) 

(3.55) 


Db  =  77b(t77||u6P  ;  Dn  =  CB(^Cj\ub\ubh  Cb  =  [ 

Jo 


dF 


dC 


dC 


(3.56) 


The  instantaneous  cross-shore  and  alongshore  velocities  are  expressed  as 


with 


and 


u{t,x,y,z)  =  U{t,x,y)  +  Ub{t,x,y)F{C) 
v{t,  X,  j/,  z)  =  V{t,  X,  y)  +  vb{t,  X,  y)F{() 

z  -I-  d{x,y) 


C  = 


h{t,x,y) 


(3.57) 

(3.58) 


(3.59) 


(3.60) 


T’  =  1  -  (3  +  0.75a)  +  aC^  for  0  <  C  <  1 

The  typical  values  of  7/^  =  0.1  and  a  =  3  are  employed  for  the  3D  computations. 
For  a  =  3, 7/2  =  0.55, 7/3  =  -0.07  and  Cb  =  15.2. 

The  time- averaged  alongshore  momentum  equation  (2.89)  is  given  by 


d 


dn 


■^dr]  1  d 


^\2 


(3.61) 


with 


S.y  =  hUV 


(3.62) 
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(3.63) 


The  time- averaged  cross-shore  energy  equation  (2.91)  is  given  by 

AjE  -|-  —(Ep)  =  —Df  —  Db 
ax 

with 

E  =  +  hU'^  +  m)  for  d>0 


E 

=  +  hU'^  +  rnj  for  d  <  0 

(3.64) 

AE 

_  E{t  =  tend)  -  E{t  -  tbegin) 

^end  ^begin 

=  TjhU  -f  +  SmU  +  m3) 

(3.65) 

Ep 

(3.66) 

Df 

—  '^hx'^b 

(3.67) 

Db 

=  CpcrCfluil^ 

(3.68) 

The  numerical  energy  dissipation  rate,  jDnumericai?  is  defined  as 

d 

T>numerical  =  -AjB  -  —{Ep)  —  Df  —  Dp  (3.69) 

which  will  be  zero  if  the  computed  results  satisfy  (3.63). 

3.2.2  Numerical  Procedures 

The  2D  numerical  model  presented  in  Section  3.1  is  extended  herein  to  in¬ 
clude  the  two  additional  equations  (3.52)  and  (3.53)  for  m  and  n,  respectively.  The 
coordinate  system,  boundary  and  initial  conditions  as  well  as  the  specified  incident 
wave  train  at  the  seaward  boundary  remain  the  same.  For  the  known  values  of 
Tj,  h  =  (t)  +  d),  C/,  m,  V  and  n  at  the  time  level  t  and  at  all  the  nodes  used  in 
the  computation,  the  values  of  these  variables  at  the  next  time  level  t*  ■=  {t  +  At), 
which  are  denoted  by  the  superscript  asterisk,  are  computed  in  sequence.  The 
variable  time  step  size  At  is  determined  for  each  time  step  using  an  approximate 
numerical  stability  criterion  for  the  explicit  finite  difference  method  adopted  in  the 
following. 


First,  along  each  of  the  cross-shore  lines,  the  values  of  /ij  =  {rjj  +  dj), 
Uj  and  mt  at  the  node  j  with  j  =  1,  2,  ■  ■  •,  s*  are  computed  from  (3.46),  (3.47) 
and  (3.52)  for  the  incident  wave  train  rn  specified  at  x  =  0.  The  integer  s*,  which 
must  be  less  than  J,  indicates  the  wet  node  next  to  the  moving  shoreline  at  the 
next  time  level  t  —  t*.  To  compute  the  cross-shore  wave  motion  along  each  line,  the 
computer  program  VBREAK  developed  for  normally  incident  waves,  as  described  in 
detail  in  the  report  by  Kobayashi  and  Johnson  (1995),  is  modified  slightly  for  the 
obliquely  incident  wave  trains.  VBREAK  solves  (3.46),  (3.47)  and  (3.52)  using  the 
MacCormack  method  (MacCormack  1969). 

The  variable  time  step  size  At  is  determined  at  the  beginning  of  each  time 

step  using  an  approximate  equation 

r*  At- 

At - ,  “  „  for  j  =  (3.70) 

max  {\Uj  \  4-  Jhj) 

in  which  Cn  is  the  Courant  number  and  the  denominator  in  (3.70)  is  the  maximum 
value  of  {\Uj\+yJhj)  at  all  the  wet  nodes  at  the  present  time  t.  The  stability  criterion 
for  the  MacCormack  method  applied  to  (3.46)  and  (3.47)  with  m  =  0  and  without 
the  terms  on  the  right  hand  side  of  (3.47)  is  Cn  <  1  {e.g.,  Anderson  et  al.  1984). 
Equation  (3.70)  is  approximate  because  the  characteristic  equations  corresponding 
to  (3.46),  (3.47)  and  (3.52)  with  m  7^  0  can  not  be  expressed  in  simple  analytical 
forms.  Moreover,  (3.70)  does  not  account  for  the  shoreline  algorithm  which  tends 
to  suffer  numerical  difficulties.  Consequently,  the  value  of  Cn  less  than  unity  is 
specified  as  input  to  adjust  At  for  successful  computation. 

Use  of  the  MacCormack  method  results  in  numerical  high-frequency  oscilla¬ 
tions  which  tends  to  appear  at  the  rear  of  a  breaking  wave,  especially  on  a  gentle 
slope  (Kobayashi  and  Johnson  1995).  The  procedure  for  smoothing  these  high- 
frequency  numerical  oscillations  presented  by  Jameson  et  al.  (1981)  and  summa¬ 
rized  by  Chaudhry  (1993)  for  hydraulic  jumps  is  modified  slightly  for  breaking  waves 
on  a  sloping  beach.  This  modified  smoothing  procedure  is  presented  in  detail  by 
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Kobayashi  and  Johnson  (1995),  who  have  stated  that  for  breaking  waves  on  gently 
sloping  beaches,  the  numerical  damping  coefficient  k  specified  as  input  to  VBREAK 
on  the  order  of  unity  is  necessary  to  damp  the  high-frequency  oscillations  whereas 
for  waves  surging  on  steep  slopes  of  coastal  structure,  the  value  of  k  on  the  order 
of  0.1  appears  to  be  sufficient.  The  computed  results  for  different  values  of  k  were 
presented  by  Johnson  et  al.  (1996). 

The  seaward  boundary  algorithm  for  and  at  x  =  0  is  based  on  the 
characteristic  equations  derived  from  (3.46)  and  (3.47)  in  a  manner  similar  to  that 
devised  by  Kobayashi  et  al.  (1987,  1989)  for  the  case  of  m  =  0.  To  obtain  the  value 
of  at  X  =  0,  (3.52)  for  m  =  C2hul  from  (3.54)  is  rewritten  in  terms  of  Uh  and 
approximated  by  an  explicit  first-order  finite  difference  based  on  the  known  values 
at  the  nodes  j  =  I  and  2.  This  is  because  m  =  0  is  a  trivial  solution  of  (3.52)  for 
X  >  0  and  t  >  0  for  the  initial  condition  m  =  0  at  t  =  0.  This  problem  could 
be  avoided  if  the  value  of  m  associated  with  the  vertical  variation  of  u  at  x  =  0 
were  known  in  the  presence  of  reflected  waves.  On  the  other  hand,  the  landward 
boundary  algorithm  dealing  with  the  moving  shoreline  is  a  minor  extension  of  that 
used  by  Kobayashi  et  al.  (1987)  for  the  case  of  m  =  0  where  the  vertical  variation 
of  u  is  assumed  to  be  small  in  the  vicinity  of  the  moving  shoreline. 

After  computing  the  values  of  rjj,  hj,  UJ  and  m’j  with  J  =  1, 2,  •  •  •,  s*  for  each 
line,  the  values  of  V^*  and  nj  at  the  node  j  with  j  =  1,  2,  s*  along  the  middle 
cross-shore  line  are  computed  using  (3.48)  and  (3.53)  which  are  expressed  in  the 
following  vector  form; 


(3.71) 


with 


R2 


(3.72) 


46 


and 


Pi  =  U*qi  +  n 

P2  =  nU*  +  m*V  +  Csh*  {ulf  Vb 

R,  =  h*^  +  fb\ul\vb 


dm* 


R2  =  fb\UbWb^b  +  fbWbl^b'^b  +  2CB<rC]\ul\ulh  -  V  — - h  n 


du* 


(3.73) 

(3.74) 

(3.75) 

(3.76) 


dx  '  dx 

where  =  hV  and  use  is  made  of  (3.50),  (3.55)  and  (3.56).  The  computed  values 
of  h*,  U*,  m*  =  C2h*{u*bY,  ul  and  ul  -  {ul  +  U*)  at  the  next  time  level  t*  along 
the  middle  cross-shore  line  i  are  used  in  (3.73)  -  (3.76).  The  term  {dr)*fdy)j  for  the 
cross-shore  line  i  at  the  node  j  is  approximated  by  a  central  finite  difference  based 
on  the  computed  values  of  rj*  at  the  node  j  at  the  two  adjacent  cross-shore  lines 
(f  —  1)  and  (z  +  1)  if  all  the  nodes  j  at  the  three  lines  (z  —  1),  z,  and  (z  +  1)  are  wet 
and  seaward  of  the  shoreline  whose  location  varies  for  the  three  cross-shore  lines. 
Otherwise,  this  term  is  set  to  be  zero. 

For  the  known  values  of  q^  =  hV,  n  =  C2hubh  from  (3.54)  and  Vb  =  {V  +  Vb) 
from  (3.51)  at  the  present  time  t  along  the  cross-shore  line,  (3.71)  is  solved  using 
the  MacCormack  method  whose  accuracy  is  second  order  in  time  and  space.  The 
predictor  and  corrector  steps  of  the  MacCormack  method  and  the  value  of  Vt  at 
the  node  j  and  at  the  next  time  level  i*  are  expressed  as 


AtR,- 


=  V,- 


Ax 


(i*i  ^j-i)  AtRj 
VJ  =  5(V,-+V,)  for  i  =  1.2,. 


(3.77) 

(3.78) 

(3.79) 


Equation  (3.77)  yields  the  temporary  predicted  values  of  q^-  and  hj.  Use  is  made 
of  Vj  =  %lh*j,  {db)j  =  nj/[C2h*j{ub)j]  and  {vb)j  =  [(v6)i  +  to  compute  Pj 
and  Rj  in  (3.78).  The  terms  dm* jdx  and  dU*/dx  in  (3.76)  are  approximated  by 
the  forward  and  backward  spatial  differences  in  (3.77)  and  (3.78),  respectively,  to 
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be  consistent  with  the  forward  and  backward  spatial  diiferences  of  dVjdx  used  in 
(3.78)  and  (3.79). 

The  moving  shoreline  is  accounted  for  in  (3.77)  by  setting  =  0  at 
j  =  (5*  +  1)  for  the  dry  node  landward  of  the  moving  shoreline.  As  for  the  seaward 
boundary  algorithm  at  the  node  j  =  1,  the  2D  model  used  the  characteristic  equa¬ 
tion  (3.42)  for  the  case  of  n  =  0.  An  alternative  and  simpler  algorithm  is  adopted 
herein  considering  the  difficulty  in  deriving  the  characteristic  equations  correspond¬ 
ing  to  (3.71).  The  imaginary  node  at  i  =  0  is  added  to  compute  Vi  using  (3.78) 
in  which  the  linearly  extrapolated  values  of  Po  =  (2Pi  -  P2),  Uq  =  {2U^  -  U2) 
and  mo  =  (2m*  -  ml)  at  the  node  i  =  0  are  employed  for  the  backward  spatial 
differences  involved  in  (3.78). 

After  q}.  and  n*  with  j  =  1,  2,  •  •  •,  s*  are  computed  using  (3.79),  =  qljh* 

is  computed  and  the  procedure  used  to  smooth  h*-,  U*  and  m^  is  applied  to  smooth 
V^*  and  n*j.  The  values  of  {vb)*j  and  (t;f,)|  are  then  computed  using  the  smoothed 
V*  and  Tij  together  with  (3.54)  and  (3.51).  Finally,  the  values  of  these  variables  for 
i  >  (-s*  +  1)  are  set  to  be  zero  for  the  dry  nodes  landward  of  the  moving  shoreline. 


Chapter  4 


COMPARISON  OF  THE  2D  MODEL  WITH 
AVAILABLE  DATA 


The  2D  numerical  model  developed  in  Section  3.1  is  compared  with  regular 
wave  laboratory  data  as  well  as  field  data.  First,  the  incident  wave  profiles  used  in 
these  comparisons  are  described  in  the  following. 

4.1  Incident  Wave  Profiles 

4.1.1  Regular  Waves 

The  incident  wave  profile  at  the  seaward  boundary  needs  to  be  specified  as 
input  to  the  numerical  model.  The  obliquely  incident  regular  wave  train  y')  at 
x'  =  0  for  the  small  angles  of  incidence  6i  in  radian  is  assumed  to  be  in  the  following 
dimensional  form: 

=  periodic  function  of  -  -^rj^  (4.1) 

in  which 

L'  =  wavelength  at  the  seaward  boundary 

C'  =  L' jT'  =  phase  velocity  at  the  seaward  boundary 

Equation  (4.1)  accounts  for  the  phase  or  time  lag  along  the  alongshore  coordinate 
y' .  The  periodic  function  in  (4.1)  is  specified  using  Stokes  second-order  or  cnoidal 


wave  theory  depending  on  the  value  of  Ursell  paranaeter  Ur  in  the  same  way  as  in 
the  existing  one-dimensional  model  (Kobayashi  et  al.  1987)  where  Ur  is  defined  as 

Ur  =  ^  (4.2) 

dt 

with 

L  =  L' ! dj  =  normalized  wavelength  at  the  seaward  boundary 
di=  d[IH'  =  normalized  water  depth  at  the  seaward  boundary 

The  periodic  function  in  (4.1)  is  specified  using  Stokes  second-order  theory  if  Ur  <  26 
and  using  cnoidal  wave  theory  if  Ur  >  26.  The  dimensional  wave  train  r)[{t\y')  is 
computed  using  (4.1)  and  normalized  using  (2.4)  and  (2.15). 


4.1.2  Irregular  Waves 

On  the  other  hand,  the  incident  random  wave  train  at  the  seaward  boundary 
is  assumed  to  be  unidirectional  and  expressible  in  the  following  dimensional  form: 

nl  (('.  y')  =  E  c;  cos  -  eXv’  +  *»  +  a*)  (4.3) 

with 


fn  = 

A/'  = 


N  = 


Oi  = 


$  = 
c*;  = 


nAf  =  frequency 
— ^  =  frequency  band  width 
computation  duration 

=  number  of  harmonics 

2Ati 

sampling  rate 

wave  number  for  bcised  on  linear  wave  theory 
representative  incident  wave  angle  {Of  <C  1) 
random  phase  angle  in  the  range  0  —  27r 
[2S'{f{JAfY^^  —  amplitude 
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S'{f^)  =  measured  incident  wave  frequency  spectrum 

A#  =  phase  shift  chosen  such  that  rj^  =  0  at  t'  =  0 
dn' 

as  well  as  >0  at  t'  =  Q  and  =  0 

The  incident  random  wave  train  t/,-  at  given  y'  is  computed  using  an  inverse 
fast  Fourier  transform  of  the  following  equation  obtained  from  (4.3): 


N 


(4.4) 


n=l 


with 


K  =  C';sin(W-4,-A«) 


The  phase  shift  A$  in  (4.3)  is  introduced  to  satisfy  the  initial  conditions  of 
Tj'i  =  0  at  f  =  0.  The  condition  of  drj'Jdt'  >  0  at  t'  =  0  is  imposed  to  produce  initial 
longshore  current  in  the  downwave  direction  for  positive  6i  for  the  coordinate  system 
shown  in  Figure  3.1.  This  procedure  reduces  the  transient  duration  associated  with 
the  slow  development  of  longshore  current  as  will  be  shown  in  Section  4.3.  To  satisfy 
the  above  two  conditions,  the  phase  shift  A$  is  computed  from  (4.3)  with  yl  =  0  at 


=  0  and  y'  =  0 


tan  A$  = 


E(r=ic;cos$. 


(4.5) 


which  has  two  roots  in  the  range  0  <  A#  <  27r.  The  condition  of  drj[ldt'  >  0  at 
t'  =  0  and  y'  =  0  requires  that 


Y,  sin($n  +  A$)  <  0  (4.6) 

n=l 

Only  one  of  the  two  roots  obtained  from  (4.5)  satisfies  (4.6). 

The  incident  random  wave  profile  rj'-  at  given  y'  is  computed  using  (4.4)  with 
(4.5)  and  (4.6).  Since  (4.5)  does  not  ensure  r}[  =  0  at  t'  =  0  for  y'  >  0,  the  computed 
values  of  rj'-  at  =  0  for  y'  >  0  are  adjusted  as  will  be  explained  in  Section  4.3. 
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4.2  Comparison  with  Regular  Wave  Data 

Visser  (1991)  conducted  eight  monochromatic  wave  experiments  in  34  m  long, 
16.6  m  wide  and  0.68  m  deep  wave  basin.  A  pumping  system  was  installed  to 
ensure  the  alongshore  uniformity  of  longshore  current.  The  beaches  were  made  of 
concrete  with  1:10  and  1:20  slopes.  Detailed  data  on  uniform  longshore  currents, 
local  wave  heights,  angles  of  wave  incidence,  wave  setup  and  runup  were  tabulated. 
The  observations  were  started  60  minutes  after  the  start  of  the  wave  generator  to 
eliminate  start-related  variations  in  flow  and  wave  fields.  Current  velocities  were 
measured  using  a  dye  method  whereas  wave  runup  was  measured  visually. 

4.2.1  Input  Parameters 

The  numerical  model  is  compared  with  four  experiments  for  which  the  sea¬ 
ward  boundary  location  can  be  taken  to  be  in  relatively  shallow  water  seaward  of 
the  breaker  line.  Table  4.1  lists  the  experiment  number  used  by  Visser  (1991)  and 
the  slope  and  incident  wave  characteristics  specified  as  input  to  the  numerical  model 
where 

tan  0'  =  uniform  slope 

d'^  =  water  depth  below  SWL  at  the  seaward  boundary 

located  at  x'  =  0 

T'  —  regular  wave  period  used  for  the  normalization 

H'  =  incident  regular  wave  height  at  x'  =  0  used  for  the  normalization 

9i  =  angle  in  degrees  of  wave  incidence  at  x'  =  0 

cr  =  ratio  between  the  cross-shore  and  vertical  length  scales 
defined  as  T'{gfH')2  in  (2.6) 

0^  =  Oi  =  reference  incident  wave  angle  in  radian  used 

for  the  normalization 

^  =  surf  similarity  parameter  (Battjes  1974)  given  by  ^ 
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Table  4.1:  Incident  waves  at  seaward  boundary  a:'  =  0. 


Expt. 

No. 

tan  9' 

d', 

(cm) 

r 

(s) 

H' 

(cm) 

Oi 

(deg) 

a 

2 

0.101 

21.1 

1.00 

9.5 

IKilfel 

0.409 

3 

0.101 

21.3 

1.00 

8.7 

14.2 

■HUM 

0.428 

4 

0.050 

18.5 

1.02 

7.9 

13.9 

11.4 

0.059 

0.227 

5 

0.050 

18.2 

1.85 

9.0 

12.9 

19.3 

0.385 

For  these  experiments,  plunging  breakers  were  observed.  The  assumptions  of 
^  1  and  <C  1  may  be  appropriate  except  for  experiment  2  with  6^  =  0.206. 
The  only  empirical  parameter  involved  in  the  numerical  model  is  the  bottom  friction 
factor  fl  in  (2.51)  where  fl  ~  0.05  has  been  used  for  predicting  wave  runup  on 
smooth  slopes  in  small-scale  experiments  {e.g.,  Kobayashi  et  al.  1989).  The  value 
of  fl  =  0.05  is  used  here  for  both  cross-shore  and  alongshore  fluid  motions. 

The  normalized  grid  sizes  Ax  and  Ay  in  Figure  3.1  need  to  be  chosen  to 
be  small  enough  to  resolve  breaking  waves  in  the  surf  and  swash  zones.  For  these 
experiments  of  alongshore  uniformity,  it  is  sufficient  to  use  the  three  cross-shore 
lines  at  y  =  0,  Ay  and  2  Ay  in  Figure  3.1.  The  value  of  Ax  is  selected  to  be  on 
the  order  of  0.01,  corresponding  to  200  grid  spacings  between  the  seaward  boundary 
and  the  still  water  shoreline.  The  value  of  Ay  is  chosen  to  be  the  same  as  Ax  to 
yield  the  same  spatial  resolution  in  the  normalized  coordinates.  A  limited  sensitivity 
analysis  with  (Ay/Ax)  ~  1, 2, 5,  and  10  has  indicated  that  the  computed  results  for 
experiment  2  remain  essentially  the  same  as  long  as  Ay  is  on  the  order  of  Ax  as  will 
be  shown  in  Figure  4.10  as  an  example. 

The  incident  wave  train  rji  as  a  function  of  time  t  at  longshore  location 
y  =  0,  Ay  and  2  Ay  is  computed  using  (4.1)  and  normalized  using  (2.4)  and  (2.5). 
The  incident  wave  train  at  y  =  0  is  computed  in  the  same  way  as  in  the  existing 
one-dimensional  model  (e.y.,  Kobayashi  et  al.  1987)  whereas  the  time  series  of  rfi  at 
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1 


y=0  - y=Ay  . y=2Ay 


Figure  4.1:  Example  of  incident  regular  wave  trains  at  j/  =  0,  Ay,  2  Ay  (alongshore 
phase  shift  is  exaggerated  for  clarity). 

y  =  Ay  and  2  Ay  are  computed  using  the  computed  time  series  of  rji{t,  0)  with  the 
alongshore  phase  shifts  depending  on  the  characteristics  of  the  incident  wave  and 
the  value  of  Ay.  As  an  example,  Figure  4.1  shows  a  Stokes  second-order  incident 
wave  train  for  experiment  2  at  location  y  —  0,  Ay  and  2  Ay  where  the  alongshore 
phase  shifts  have  been  exaggerated  for  clarity. 

4.2.2  Computed  Time  Series 

The  detailed  computed  results  for  experiment  2  are  presented  as  an  example. 
The  temporal  and  cross-shore  variations  of  the  free  surface  elevation  r],  the  depth- 
averaged  cross-shore  velocity  U,  and  the  depth-averaged  alongshore  velocity  V  are 
stored  along  the  center  line  &t  y  =  Ay. 

The  temporal  variations  of  r],  U  and  V  for  the  duration  0  <  t  <  300  at 
x  =  0  (at  the  seaward  boundary),  x  =  0.509  (immediately  seaward  of  the  breaker 
line),  X  =  0.770  (in  the  outer  surf  zone),  x  =  1.550  (in  the  inner  surf  zone),  and 
X  =  2.265  (in  the  swash  zone)  are  shown  in  Figures  4.2,  4.3  and  4.4,  respectively.  The 
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cross-shore  fluid  motion  represented  by  Tj  and  U  computed  using  (3.12)  and  (3.13) 
becomes  periodic  fairly  quickly  for  t  >  20  as  has  been  the  case  with  the  previous  one¬ 
dimensional  computations  for  beaches  {e.g.,  Kobayashi  et  al.  1989).  The  alongshore 
fluid  motion  represented  by  V  computed  using  (3.34)  becomes  periodic  very  slowly 
especially  in  the  vicinity  of  the  breaker  line.  The  very  slow  response  of  the  alongshore 
fluid  motion  is  qualitatively  consistent  with  the  analytical  result  of  Ryrie  (1983)  for 
the  periodic  solution  development  as  well  as  the  experimental  procedure  adopted  by 
Visser  (1991)  who  made  measurements  one  hour  after  the  start  of  the  wave  maker. 
In  light  of  Figure  4.4,  the  time  averaging  denoted  by  the  overbar  in  the  following  is 
performed  for  the  duration  200  <  t  <  300. 

The  computed  cross-shore  variations  of  rj,  j/rms,  U,  Urms,  V  and  I4ms  for 
experiment  2  are  shown  in  Figure  4.5  where  the  root-mean-square  values  representing 
the  magnitude  of  the  oscillatory  components  are  defined  as 

=  ;  u^„„  =  {u-uy  ;  KL  =  (v'-V)'  (4,7) 

For  example,  if  ?;  =  0.5cos(27rt),  ^  =  0  and  T/rms  =  l/\/8  =  0.35.  The  normalized 
uniform  slope  is  indicated  by  the  dashed-dotted  straight  line  in  the  top  panel  in 
Figure  4.5.  The  upper  limit  of  the  wave  setup  fj  is  the  maximum  runup  elevation 
on  the  slope  above  SWL  because  h  >  0  in  the  region  wetted  by  water.  The  increase 
of  Tf  and  the  decrease  of  rj^ms  in  the  surf  and  swash  zones  are  approximately  linear. 
On  the  other  hand,  Urms  decreases  slowly  in  the  surf  zone  and  rapidly  in  the  swash 
zone.  U  is  negative  and  represents  the  cross-shore  return  current  as  explained  by 
Kobayashi  et  al.  (1989).  The  longshore  current  V  is  dominant  in  the  surf  zone 
and  the  oscillatory  component  decreases  approximately  linearly  in  the  surf  and 
s Welsh  zones. 

Figure  4.6  shows  the  computed  cross-shore  variations  of  the  normalized  quan¬ 
tities  involved  in  the  time-averaged  cross-shore  wave  energy  equation  (3.7)  and  the 
time-averaged  alongshore  momentum  equation  (3.5).  In  Figure  4.6,  E  =  specific 
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energy  per  unit  horizontal  area;  Ep  =  energy  flux  per  unit  width;  Dj  =  energy 
dissipation  rate  due  to  bottom  friction  per  unit  horizontal  area;  Dp  =  energy  dissi¬ 
pation  rate  due  to  wave  breaking  per  unit  horizontal  area;  Sxy  =  alongshore  radiation 
stress  given  by  (3.6);  and  fb\U\V  =  alongshore  bottom  shear  stress. 

Figure  4.6  indicates  that  the  energy  dissipation  due  to  wave  breaking  is  dom¬ 
inant  and  does  not  occur  suddenly  in  this  numerical  model  which  does  not  account 
for  wave  breaking  explicitly  (Kobayashi  and  Wurjanto  1992).  For  these  experiments 
of  alongshore  uniformity,  the  computed  alongshore  gradients  of  the  mean  and  vari¬ 
ance  of  7/  are  negligible  and  the  time-averaged  alongshore  momentum  equation  (3.5) 
reduces  to  dSxyjdx  =  —fb\U\V.  Figure  4.6  also  shows  that  Sxy  decreases  monoton- 
ically  in  the  surf  and  swash  zones.  The  computed  cross-shore  variations  of  dSxyldx 
and  -fb\U\V  are  essentially  the  same  where  fb\U\V  is  plotted  to  distinguish  the  two 
curves. 

4.2.3  Comparison  with  Measurements 

Figure  4.7  compares  the  measured  and  computed  cross-shore  variations  of  the 
normalized  local  wave  height  H  for  each  of  the  four  experiments  listed  in  Table  4.1. 
The  agreement  is  very  good  in  view  of  no  adjustable  parameter  included  in  this 
numerical  model  to  initiate  wave  breaking.  However,  it  should  be  stated  that  this 
numerical  model  can  not  predict  wave  shoaling  without  wave  breaking  over  a  long 
distance  (Kobayashi  et  al.  1989). 

Figure  4.8  compares  the  measured  and  computed  cross-shore  variations  of 
the  normalized  wave  setup  rj  together  with  the  normalized  uniform  slope  indicated 
by  the  dashed  line  for  each  of  the  four  experiments.  The  agreement  is  good  in 
the  swash  zone  but  the  computed  mean  water  level  rises  too  rapidly  landward  of 
the  breaker  line  as  was  the  case  with  the  previous  comparison  by  Kobayashi  et  al. 
(1989).  The  numerical  model  does  not  predict  the  transition  zone  of  constant  wave 


set-down  whose  effects  on  surf  zone  hydrodynamics  were  reviewed  and  elaborated 
by  Nairn  et  al.  (1990). 

Table  4.2  shows  the  comparisons  of  the  measured  and  computed  maximum 
setup  and  runup  for  the  four  experiments.  The  computed  maximum  setup  and  runup 
correspond  to  the  mean  and  maximum  shoreline  elevations,  respectively,  measured 
by  hypothetical  wires  placed  parallel  to  and  above  the  uniform  slope  at  elevations 
of  1,  5  and  10  mm,  whereas  the  actual  measurements  were  made  visually. 

The  computed  maximum  setup  and  runup  are  not  very  sensitive  to  the  wire 
elevations  and  in  fair  agreement  with  the  measured  values  except  that  the  numerical 
model  with  the  bottom  friction  factor  fl  =  0.05  slightly  underpredicts  the  visually 
measured  runup.  It  is  also  noted  that  the  swash  oscillations  in  the  regular  wave 
experiments  are  very  narrow  in  comparison  to  swash  oscillations  on  natural  beaches 
that  tend  to  be  dominated  by  low-frequency  motions  (Guza  and  Thornton  1982, 
Holman  and  Sallenger  1985). 


Table  4.2:  Measured  and  computed  maximum  setup  and  runup. 


Expt. 

No. 

Maximum  Setu] 

3 

Runup 

Computed 

Meas. 

Computed 

Meas. 

Imm 

5mm 

10mm 

1mm 

5mm 

10mm 

2 

0.34 

0.29 

0.29 

0.38 

0.35 

0.35 

0.43 

3 

0.34 

0.28 

MM 

0.31 

0.38 

0.35 

0.35 

0.47 

4 

0.19 

0.17 

0.20 

0.20 

0.20 

0.21 

0.24 

5 

0.28 

0.23 

0.27 

0.31 

0.29 

0.29 

0.34 

Figure  4.9  compares  the  measured  and  computed  cross-shore  variations  of  the 
longshore  current  V  for  the  four  experiments.  The  numerical  model  with  fl  =  0.05 
predicts  the  magnitude  of  V  but  can  not  predict  the  shape  of  V  probably  because 
the  numerical  model  based  on  (3.3)  does  not  include  lateral  mixing  (dispersion) 
and  it  can  not  predict  the  transition  zone  as  shown  in  Figure  4.8.  Comparing 
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1.0 


Experiment  2 


Computed  temporal  variations  of  free  surface  elevation  t]  at  x=0,  0.509, 
0.770,  1.550  and  2.265. 
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Experiment  2 


Computed  temporal  variations  of  depth- averaged  alongshore  velocity 
V  aXx  =  0,0.509,0.770,1.550  and  2.265. 
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Experiment  2 


Figure  4.6:  Time-averaged  cross-shore  wave  energy  and  alongshore  momentum 
balances. 


63 


X 


Figure  4.9:  Measured  and  computed  longshore  current  V  for  four  experiments. 


the  similar  agreement  for  experiments  2  and  3  whose  incident  wave  conditions  are 
listed  in  Table  4.1,  it  may  be  concluded  that  0^  =  0.206  may  still  be  regarded  to 
be  much  less  than  unity.  Visser  (1984)  and  Nairn  et  al.  (1990)  showed  it  would 
be  necessary  to  delay  the  initiation  of  the  influence  of  energy  dissipation  on  the 
generation  of  longshore  currents  until  the  landward  limit  of  the  transition  zone. 
These  shortcomings  of  the  numerical  model  may  be  serious  for  longshore  currents 
generated  by  regular  waves  but  are  much  less  apparent  for  irregular  waves  due  to 
irregular  wave  breaking  and  generation  of  low-frequency  motions  as  can  be  seen  in 
the  comparison  with  the  field  data  of  Thornton  and  Guza  (1986)  in  the  next  section. 

Finally,  Figure  4.10  indicates  that  the  computed  longshore  current  is  not 
sensitive  to  the  values  of  (Ay/ Ax)  used  in  the  computation  as  long  as  Ay  is  on  the 
order  of  Ax  where  Ay  ~  Ax  has  been  used  unless  stated  otherwise. 

4.3  Comparison  with  Field  Data 

The  2D  numerical  model  is  also  compared  with  field  data  to  verify  whether 
the  numerical  model  is  applicable  to  natural  beaches  as  well.  The  numerical  model 
is  compared  herein  with  the  field  experiments  conducted  by  Thornton  and  Guza 
(1986,  1989)  at  Leadbetter  beach  in  February  1980.  The  mean  nearshore  slope 

varied  between  0.03  and  0.06  during  the  experiment  but  no  major  bar  was 
apparent.  The  incident  waves  were  limited  to  a  narrow  window  of  approach  because 
of  the  protection  from  a  cape  and  islands.  Comparison  is  made  with  the  data  of 
February  5  and  6  that  included  the  bottom  profiles,  the  incident  wave  frequency 
spectrum  and  representative  angle  of  incidence  in  3  m  depth,  and  the  cross-shore 
variations  of  the  measured  root-mean-square  wave  height  and  longshore  current. 
The  wave  conditions  and  bottom  contours  near  the  shoreline  could  be  reasonably 
assumed  to  be  uniform  in  the  alongshore  direction. 


Figure  4.10:  Sensitivity  analysis  for  experiment  2  with  (a)  Ay 
(b)  Ay  ~  2Aa:,  (c)  Ay  ~  5Aa:,  (d)  Ay  ~  lOAx. 


4.3.1  Input  Parameters 

The  measured  incident  wave  frequency  spectra  in  the  3  m  water  depth  on 
February  5  and  6  presented  by  Thornton  and  Guza  (1986)  included  only  the  sea-swell 
band  of  frequencies  (0.05-0.3  Hz).  As  a  result,  the  incident  wave  train  computed 
using  (4.4)  does  not  include  the  incident  low-frequency  wave  components.  The  spec¬ 
tral  peak  period  the  root-mean-square  wave  height  and  the  representative 

incident  wave  angle  9i  at  the  seaward  boundary  in  the  depth  =  3  m  are  summa¬ 
rized  in  Table  4.3.  The  assumption  of  random  phases  in  (4.4)  is  not  really  valid  in 
shallow  water  due  to  nonlinear  phase  coupling  {e.g.,  Elgar  and  Guza  1986).  The 
computed  results  presented  in  the  following  will  need  to  be  interpreted  in  light  of 
these  limitations  of  (4.4). 

The  computed  results  presented  in  the  following  are  based  on  the  normal¬ 
ization  using  H'  =  T'  =  and  Oc  =  Oi  in  radians  at  the  3  m  depth.  The 

parameter  a  defined  in  (2.6)  is  hence  given  by 


1/2 


H' 

K  rms/ 


The  values  of  a  and  Oc  listed  in  Table  4.3  suggest  that  the  assumptions  of  shallow 
water  waves  with  small  angles  of  incidence  are  appropriate  for  these  data  sets. 


Table  4.3:  Incident  wave  at  seaward  boundary  d'  =  3  m. 


Date 

d' 

(m) 

t; 

(s) 

IJf 

^rms 

(m) 

0i 

(deg) 

<7 

0^c 

Feb  5 

■niK« 

3.0 

12.8 

7.8 

57 

0.0196 

0.80 

Feb  6 

3.0 

11.1 

IS 

7.5 

66 

0.0169 

0.87 

The  number  of  cross-shore  lines  in  the  finite  difference  grid  shown  in  Fig¬ 
ure  3.1  is  taken  as  7  =  3  since  the  beach  and  incident  wave  conditions  are  assumed 
to  be  uniform  in  the  longshore  direction.  The  bottom  profiles  along  the  three  cross¬ 
shore  lines  are  assumed  to  be  the  same  and  taken  as  the  measured  bottom  profiles 
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on  February  5  and  6.  The  mean  bottom  slope,  tan^^,  between  the  shoreline  and 
the  mean  breaker  line  was  0.035  and  0.033  on  February  5  and  6,  respectively,  as 
listed  in  Table  4.3.  The  actual  bottom  profiles  shown  in  Figure  4.11  are  used  for 
the  following  computations.  Plunging  breakers  were  most  often  observed.  The  surf 
similarity  parameter  defined  as  ^  =  atan  0^/(27r)^/^  was  ^  =  0.80  and  0.87  on 
February  5  and  6. 

The  dimensional  cross-shore  grid  spacing  is  taken  as  Ax'  =  0.180  m  on  Febru¬ 
ary  5  with  =  0.49  m  and  Ax'  =  0.135  m  on  February  6  with  =  0.28  m  in 
order  to  resolve  steep  wave  fronts.  The  normalized  grid  spacing  Ax  is  about  0.007 
on  both  days.  The  number  of  grid  spacings  between  the  3  m  depth  and  the  still 
water  shoreline  is  400  and  520  on  February  5  and  6,  respectively.  The  dimensional 
alongshore  grid  spacing  is  taken  as  Ay'  =  1.32  and  1.03  m  on  February  5  and  6, 
respectively,  so  that  the  normalized  grid  spacings  Ax  and  Ay  are  approximately  the 
same. 

The  normalized  computation  duration,  tra&x  =  ^max/^p?  taken  as  fmax  =  500. 
Correspondingly,  =  107  and  93  min  on  February  5  and  6,  respectively.  The 


computation  duration  exceeds  the  sampling  duration  of  68  min  employed  by 
Thornton  and  Guza  (1986)  to  account  for  the  transient  duration  of  the  computation 
initiated  from  no  wave  action  at  t  =  0. 

The  incident  wave  train  computed  at  y'  =  0,  Ay'  and  2A?/'  using 

(4.4)  is  normalized  as  rji  =  with  t  =  t' /T^.  The  sampling  rate  At'  associated 

with  N  in  (4.4)  was  taken  to  be  0.5  sec.  The  normalized  sampling  rate  Ats  must  be 
small  enough  to  resolve  the  temporal  variation  of  rji  but  is  normally  much  larger  than 
the  finite  difference  time  step  At  required  by  the  numerical  stability  criterion.  The 
time  step  At  used  for  the  actual  computation  is  about  0.0002  and  much  smaller  than 
the  sampling  rate  to  minimize  spurious  shoreline  oscillations  as  will  be  explained 
later  in  the  discussion  of  the  computation  results  for  alongshore  velocity.  A  simple 
linear  interpolation  of  r}i  sampled  at  the  rate  At^  is  performed  to  find  the  value  of  rji 
at  the  time  level  t*  =  {t  +  At)  during  the  time-marching  computation.  The  values 
of  Tji  at  t  =  0  is  set  to  zero  for  y'  >  0  since  (4.5)  does  not  ensure  =  0  at  t'  =  0 
and  y'  >  0.  The  effect  of  these  adjustment  to  the  specified  spectrum  is  extremely 
small.  Figure  4.12  shows  the  normalized  incident  wave  trains  rji  aX  y  =  0,  Ay  and 
2 Ay  which  satisfy  the  initial  condition  r/j  =  0  at  time  t  =  0,  and  drjifdt  >  0  for 
small  time  t. 

The  only  empirical  parameter  involved  in  the  numerical  model  is  the  bottom 
friction  factor  fl  in  (3.4).  Use  is  made  of  fl  =  0.015  which  was  the  value  calibrated 
by  Raubenheimer  et  al.  (1995)  using  runup  measurements  on  the  fine-grained, 
gently  sloping  Scripps  beach.  The  Leadbetter  beach  is  composed  of  fine  to  medium 
size  sand.  Thornton  and  Guza  (1986)  calibrated  the  friction  factor  Cj  for  their 
longshore  current  model  using  the  measured  longshore  currents  on  the  Leadbetter 
beach.  The  two  friction  factors  are  related  by  (7/  =  0.5/^',  assuming  that  the  bottom 
friction  factor  for  the  present  time-dependent  model  and  their  time-averaged  model 
for  longshore  current  are  the  same.  Their  calibrated  value  of  (7/  ~  0.006  turns 


out  to  be  very  consistent  with  fl  ~  0.015  based  on  the  runup  measurements  by 
Raubenheimer  et  al.  (1995). 


Figure  4.12:  Initial  portions  of  adjusted  incident  random  wave  trains  rji  at  y  =  0 
(line  1),  Ay  (line  2),  2 Ay  (line  3). 

Figure  4.13  shows  the  measured  and  normalized  incident  wave  frequency 
spectrum  Si{f)  for  February  5  and  6  corresponding  to  the  incident  wave  train  y) 
at  X  =  0  and  at  y  =  0,  Ay  and  2Ay  specified  as  input  to  the  numerical  model, 
where  the  dimensional  frequency  f  is  normalized  as  /  =  f'T^  and  the  peak  of  Si{f) 
is  located  at  /  =  1.  The  spectra  Si{f)  at  y  =  0,  Ay  and  2 Ay  are  identical  because 
of  the  assumption  of  alongshore  uniformity  of  the  incident  wave  conditions.  The 
specified  spectra  Si{f)  do  not  include  any  low-frequency  component  where  Thornton 
and  Guza  (1986)  used  f  =  0.05  Hz  as  a  cutoff  frequency.  The  normalized  reflected 
wave  trains  rjr{t.,y)  at  x  =  0  and  at  y  =  0,  Ay  and  2Ay  are  computed  when  h  and 
U  are  computed  using  (3.12)  and  (3.13)  (Kobayashi  et  al.  1989).  The  corresponding 
reflected  wave  spectra  are  computed  to  be  essentially  identical  and  are  also  shown  in 
Figure  4.13  which  present  the  smoothed  spectra  Sr(f)  with  64  degrees  of  freedom. 


The  computed  reflected  wave  spectra  5'r(/),  especially  for  February  5,  are  domi¬ 
nated  by  the  low-frequency  wave  components.  The  cross  spectra  among  r}r{t,y)  at 
y  =  0,  Ay  and  2Ay  are  also  computed.  The  computed  coherence  squared  is  nearly 
unity,  and  the  computed  phase  difference  is  almost  zero  because  of  the  very  small 
value  of  Ay  ~  0.007  used  in  the  computation.  The  reflected  wave  direction  may 
become  discernible  if  the  number  of  cross-shore  lines  and  the  value  of  Ay  are  in¬ 
creased. 

4.3.2  Computed  Time  Series 

The  computed  temporal  and  cross-shore  variations  of  t],  U  and  V  aX  y  =  Ay 
only  are  used  because  of  the  assumed  alongshore  uniformity.  Figures  4.14,  4.15,  and 
4.16  show  the  computed  temporal  variations  of  the  free  surface  elevation  y  above 
SWL,  the  depth-averaged  cross-shore  velocity  U  and  the  depth-averaged  alongshore 
velocity  V,  respectively,  at  x  =  0  (in  the  3  m  depth),  a:  =  1.28  (in  the  breaker  zone), 
X  =  1.92  (in  the  inner  surf  zone),  x  =  2.56  (near  the  still  water  shoreline),  and 
X  =  2.81  (above  the  still  water  shoreline)  for  February  5  while  Figures  4.19,  4.20, 
and  4.21  show  the  computed  temporal  variations  of  rj,  U,  and  F  at  x  =  0  (in  the 
3  m  depth),  x  =  1.98  (in  the  breaker  zone),  x  =  3.08  (in  the  inner  surf  zone), 
X  =  3.82  (near  the  still  water  shoreline),  and  x  =  3.93  (above  the  still  water  shoreline) 
for  February  6.  These  figures  are  explained  after  some  of  the  details  associated  with 
the  numerical  method  are  discussed  below. 

Figure  4.17  shows  that  the  computed  alongshore  fluid  motions  without  the 
numerical  damping  term  Dj  given  by  (3.38),  where  the  numerical  damping  coeffi¬ 
cient  Ci  is  set  to  be  zero,  are  practically  the  same  as  those  computed  using  =  1 
except  for  the  slight  modifications  of  the  high-frequency  components.  Computa¬ 
tion  is  also  made  by  reversing  the  spatial  differencing  in  the  predictor  and  corrector 
equations  (3.35)  and  (3.36).  Figure  4.18  shows  that  the  computed  results  are  in¬ 
distinguishable,  implying  that  the  MacCormack  method  adopted  in  Section  3.1.2 
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should  be  applicable  to  wave  fronts  propagating  seaward  as  well  as  landward. 

Figures  4.14  and  4.19  and  Figures  4.15  and  4.20  indicate  that  the  transient 
duration  of  tj  and  U  computed  using  (3.12)  and  (3.13)  with  no  wave  action  at  f  =  0 
is  short  as  has  been  the  case  with  the  computed  results  for  regular  waves  shown 
in  Figures  4.2  and  4.3.  On  the  other  hand,  Figures  4.16  and  4.21  show  that  the 
transient  duration  of  V  computed  using  (3.34)  is  long  near  the  breaker  zone  and 
inner  surf  zone  where  the  longshore  current  develops  very  slowly.  The  slow  response 
of  the  alongshore  fluid  motion  for  irregular  waves  is  qualitatively  similar  to  that 
in  Figure  4.4  for  regular  waves.  The  oscillatory  components  of  7/  in  Figures  4.14 
decrease  landward  where  the  lower  limit  of  r;  at  a;  =  2.56  and  2.81  is  limited  by  the 
bottom  elevation.  On  the  other  hand,  the  oscillatory  components  of  U  and  V  in 
Figures  4.15  and  4.16  increase  landward  and  become  the  maximum  near  the  still 
water  shoreline.  The  comparison  of  the  oscillatory  components  of  77,  f/  and  V  at 
X  =  2.81  above  the  still  water  shoreline  indicates  that  the  fluid  velocities  U  and  V 
oscillate  much  more  than  the  water  depth  h  in  the  swash  zone  where  h  equals  the 
free  surface  elevation  tj  above  the  bottom  elevation.  The  results  for  February  6  are 
also  similar  as  depicted  in  Figures  4.19,  4.20,  and  4.21. 

For  the  computation  for  February  6,  it  is  found  that  the  value  of  At  smaller 
than  that  required  by  the  numerical  stability  criterion  is  necessary  to  minimize  spu¬ 
rious  shoreline  oscillations.  The  computed  temporal  variations  of  alongshore  velocity 
V  using  At  ~  0.0007  plotted  in  Figures  4.22  and  4.23  show  sudden  jumps  which 
originate  at  the  shoreline  and  propagate  in  the  seaward  direction.  The  shoreline  os¬ 
cillations  along  the  three  cross-shore  lines  depicted  in  Figure  4.24  confirms  that  the 
sudden  jumps  are  caused  by  spurious  large  shoreline  differences  lasting  for  a  short 
duration.  The  use  of  At  equal  to  one  half  of  0.0007  reduces  the  jumps  but  does 
not  eliminate  the  problem  completely  as  can  be  seen  in  Figure  4.25.  Therefore,  the 
computation  is  repeated  with  At  ~  0.0002  and  the  numerical  algorithm  dealing  with 


the  moving  shoreline  is  improved  somewhat  to  minimize  spurious  shoreline  oscilla¬ 
tions.  This  problem  was  not  apparent  in  the  previous  one- dimensional  cross-shore 
computations  because  their  effects  on  the  cross-shore  fluid  motion  were  limited  to 
the  swash  zone  and  lasted  only  for  a  short  duration.  The  computed  results  for  rj  and 
U  using  these  three  different  values  of  At  remain  essentially  the  same.  As  a  result, 
the  computations  for  both  February  5  and  6  presented  in  this  chapter  are  based  on 
At  ~  0.0002  with  no  apparent  unrealistic  shoreline  oscillation. 

4.3.3  Computed  Spectra 

The  computed  alongshore  velocities  V  in  the  region  of  the  slow  longshore  cur¬ 
rent  development  in  Figures  4.16  and  4.21  exhibit  very  slow  oscillations  which  will 
be  elaborated  in  the  following.  The  computed  temporal  variations  of  ?/,  U  and  V  for 
the  duration  200  <t<  500  are  used  in  the  following  spectral  and  statistical  calcula¬ 
tions  to  account  for  the  slow  development  of  the  alongshore  velocity  V.  Admittedly, 
the  time-dependent  model  is  not  efficient  computationally  for  the  alongshore  fluid 
motion  because  of  its  slow  response  and  dominant  longshore  current. 

Figures  4.26,  4.27,  and  4.28  show  the  smoothed  frequency  spectra  of  77,  t/,  and 
V  with  64  degrees  of  freedom  for  February  5  that  correspond  to  the  time  series  of  7/,  17, 
and  V  for  200  <  t  <  500  shown  in  Figures  4.14,  4.15,  and  4.16,  respectively.  The 
noticeable  low-frequency  (/  <  0.64  corresponding  to  /'  =  0.05  Hz)  wave  components 
near  the  shoreline  are  expected  from  the  previous  field  measurements  {e.g.,  Guza  and 
Thornton,  1982;  Holman  and  Sallenger,  1985)  and  will  be  discussed  further  at  the 
end  of  this  section.  Figure  4.28  clearly  shows  the  very  low  frequency  component  of  V 
at  jc  =  1.28  (in  the  breaker  zone)  and  1.92  (in  the  inner  surf  zone)  corresponding  to 
the  very  slow  oscillations  of  V  in  Figure  4.16.  The  corresponding  very  low  frequency 
components  of  rj  and  U  at  x  =  1.28  and  1.92  are  absent  or  extremely  small  in 
Figures  4.26  and  4.27. 
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Figure  4.16: 


Computed  temporal  variations  of  depth-averaged  alongshore  velocity 
1/  at  X  =  0,  1.28,  1.92,  2.56  and  2.81  for  February  5. 
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February  5 


Figure  4.18:  Comparison  of  computed  temporal  variations  of  depth- averaged 
alongshore  velocity  at  x  =  0,  1.28,  1.92,  2.56  and  2.81  for  Febru¬ 
ary  5  using  original  scheme  and  by  reversing  the  spatial  dilferenc- 
ing  in  predictor  and  corrector  equations  of  MacCormack  method  for 
t=480-500. 
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Figure  4.21:  Computed  temporal  variations  of  depth-averaged  alongshore  velocity 
V  at  X  =  0,  1.98,  3.08,  3.82  and  3.93  for  February  6. 
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Figure  4.22:  Computed  temporal  variations  with  sudden  jumps  of  depth-averaged 
alongshore  velocity  V  at  x  =  0,  1.98,  3.08,  3.82  and  3.93  for  Febru¬ 
ary  6  with  At  ci  0.0007. 
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Figure  4.23:  Computed  temporal  variations  with  sudden  jumps  of  depth-averaged 
alongshore  velocity  F  at  x  =  0,  1.98,  3.08,  3.82  and  3.93  for  Febru¬ 
ary  6  with  At  ~  0.0007  plotted  from  t  =  50  to  t  =  100. 
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Figure  4.25;  Computed  temporal  variations  with  reduced  jumps  of  depth-averaged 
alongshore  velocity  at  x  =  0,  1.98,  3.08,  3.82  and  3.93  for  Febru¬ 
ary  6  with  At  ~  0.00035. 


87 


To  infer  the  likely  origin  of  the  very  low  frequency  component  of  V,  the 
corresponding  spectra  of  dT]ldy,qiU,  and  \U\V  are  plotted  in  Figures  4.29,  4.30,  and 
4.31,  respectively.  The  alongshore  free  surface  gradient  dTjjdy  in  the  approximate 
alongshore  momentum  equation  (3.34)  or  (3.42)  drives  the  alongshore  wave  motion 
in  this  simplified  numerical  model.  The  alongshore  momentum  flux,  hUV  =  qeJJ  in 
(3.34)  affects  the  alongshore  fluid  motion  through  the  nonlinear  interaction  of  h,  17, 
and  V,  where  h  is  the  sum  of  rj  and  the  still  water  depth.  The  alongshore  bottom 
stress  fb\U\V  in  (3.34)  reduces  the  alongshore  fluid  motion,  where  the  normalized 
bottom  friction  factor  defined  in  (3.4)  is  fi  =  0.43  on  February  5.  Figure  4.29 
indicates  the  considerably  very  low  frequency  component  of  dqjdy  at  x  =  1.28 
and  1.92  which  does  not  vary  much  with  the  normalized  frequency  /.  Comparing 
Figures  4.27  and  4.30,  the  spectral  shapes  of  U  and  q^U  at  x  =  1.28, 1.92,  and  2.56  are 
similar  in  the  region  where  the  mean  is  a  significant  portion  of  q^.  The  temporal 
variations  of  qe  are  similar  to  those  of  V  shown  in  Figure  4.16.  The  cross-shore 
variation  of  will  be  presented  later  in  relation  to  that  of  the  longshore  current  V. 
Comparison  of  Figures  4.28  and  4.31  indicates  that  the  spectral  shapes  of  V  and 
\U\V  are  similar  except  near  the  shoreline  where  the  spectra  of  \U\V  contain  more 
higher-frequency  (/  >  1.5)  components,  probably  due  to  the  nonlinear  interaction  of 
\U\  and  y. 

To  examine  whether  the  very  low  frequency  component  of  dqjdy  in  Fig¬ 
ure  4.29  is  capable  of  generating  the  very  low  frequency  component  of  V  shown  in 
Figure  4.28,  use  is  made  of  the  linearized  alongshore  momentum  equation  corre¬ 
sponding  to  (3.42) 


(4.8) 


dV  dq 
dt  ~  dy 

The  temporal  variation  of  V  for  200  <  t  <  500  is  calculated  by  integrating  (4.8) 
for  the  computed  temporal  variation  of  dqfdy  at  given  x  corresponding  to  each 
spectrum  shown  in  Figure  4.29.  The  initial  value  of  y  at  t  =  200  is  taken  as  the 


value  of  y  at  t  =  200  computed  using  (3.34). 

Figure  4.32  compares  the  spectrum  of  the  temporal  variation  V  at  given  x 
calculated  using  the  linear  equation  (4.8)  with  the  corresponding  spectrum  of  V 
computed  using  the  nonlinear  equation  (3.34)  and  shown  in  Figure  4.28.  It  is  noted 
that  the  linear  equation  (4.8)  amplifies  the  very  low  frequency  component  of  V 
because  (4.8)  yields  |y|  =  \A\/{27rf),  where  V  and  A  are  the  complex  amplitude  of 
V  and  d’q/dy  for  given  /,  respectively.  Figure  4.32  shows  that  the  two  spectra  of 
y  at  re  =  0  are  essentially  identical  except  for  the  very  low  frequency  range.  This 
indicates  that  the  nonlinear  terms  in  (3.34)  are  negligible  at  re  =  0  as  may  also  be 
inferred  from  the  spectra  of  qiU  and  |I7|y  at  re  =  0  shown  in  Figures  4.30  and  4.31. 
Figure  4.32  also  indicates  that  the  nonlinear  terms  in  (3.34)  reduce  the  spectra  of 
y  at  re  =  1.28  and  1.92,  especially  for  the  low-frequency  range.  On  the  other  hand, 
the  term  drjidy  in  (3.34)  and  (3.42)  is  negligible  relative  to  the  nonlinear  terms  near 
the  shoreline  at  rc  =  2.56  and  2.81  except  for  the  very  low  frequency  range.  This  is 
consistent  with  the  significant  decrease  of  the  spectra  of  dy/dy  near  the  shoreline 
as  shown  in  Figure  4.29. 

Figure  4.32  suggests  that  the  very  low  frequency  component  of  drjjdy  is  large 
enough  to  cause  the  corresponding  very  slow  oscillations  of  V.  However,  the  origin 
of  this  very  low  frequency  component  of  dy / dy  is  not  certain.  It  has  been  mentioned 
that  the  earlier  computations  shown  in  Figures  4.22-4.25  exhibited  infrequent  sud¬ 
den  jumps  in  the  computed  temporal  variations  of  V  caused  by  the  infrequent  large 
differences  of  the  shoreline  location  along  the  cross-shore  lines  1  and  3.  The  cen¬ 
tral  finite  difference  approximation  of  dyjdy  along  the  cross-shore  line  2  using  the 
computed  values  of  y  along  the  two  adjacent  cross-shore  lines  is  sensitive  to  the 
spurious  large  shoreline  difference  lasting  for  a  short  duration.  The  final  results  pre¬ 
sented  herein  are  computed  using  the  improved  numerical  algorithm  dealing  with 
the  moving  shoreline  and  using  a  smaller  time  step  At  as  explained  in  relation  to 
Figures  4.22-4.25.  The  computed  shoreline  oscillations  along  the  three  cross-shore 


lines  do  not  show  any  unrealistic  shoreline  location  dilference.  The  very  low  fre¬ 
quency  component  of  the  spectrum  of  drj/dy  &t  x  =  2.81  shown  in  Figure  4.29  is 
indeed  extremely  small. 

The  very  low  frequency  component  of  drildy  seems  to  be  related  to  the  break¬ 
ing  of  obliquely  incident  random  waves  because  the  very  low  frequency  components 
of  drjldy  and  V  tend  to  occur  in  the  vicinity  of  the  breaker  zone  as  shown  in  Fig¬ 
ures  4.16,  4.28,  and  4.29.  No  very  slow  oscillations  were  found  in  the  computed 
temporal  variations  of  V  for  regular  waves  as  shown  in  Figure  4.4.  Cross-spectral 
analyses  of  V  and  drjldy  do  not  indicate  any  cross-shore  wave  pattern  in  these  very 
low  frequency  components.  Since  the  present  numerical  model  does  not  simulate 
wave  breaking  explicitly,  the  origin  of  the  very  low  frequency  component  of  drj/dy 
can  still  be  numerical. 

4.3.4  Comparison  with  Measurements 

Figure  4.33  shows  the  comparisons  between  the  measured  and  computed 
cross-shore  variations  of  the  local  root-mean-square  wave  height,  JTrms  =  H\ 

on  February  5  and  6.  Thornton  and  Guza  (1986)  obtained  H!^  using  a  zero  up- 
crossing  method  from  68  minutes  free  surface  elevation  records  after  band-pass  fil¬ 
tering  between  0.05  and  0.5  Hz.  Accordingly,  the  computed  temporal  variation  of  rj 
for  200  <  i  <  500  is  analyzed  in  the  same  way  to  obtain  i/rms  based  on  the  zero  up- 
crossing  method.  In  addition,  the  computed  spectral  estimate  of  Hrms  =  (8mo)^/^ 
is  also  presented  in  Figure  4.33  where  is  the  zero  moment  of  the  frequency  spec¬ 
trum  of  rj  after  the  band-pass  filtering.  The  agreement  between  the  data  and  the 
present  time-dependent  model  is  no  better  than  that  between  the  data  and  the  time- 
averaged  model  of  Thornton  and  Guza  (1983)  presented  with  their  data  except  that 
the  only  empirical  parameter  involved  in  the  time-dependent  model  is  the  bottom 
friction  factor  whose  effect  on  the  free  surface  elevation  rj  is  essentially  limited  in 
the  swash  zone  {e.g.,  Kobayashi  and  Wurjanto  1992). 
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Figure  4.26:  Computed  frequency  spectra  of  free  surface  elevation  rj  at 
X  =  0,1.28,1.92,2.56,  and  2.81  on  February  5. 


Figure  4.28:  Computed  frequency  spectra  of  alongshore  velocity  V  at 
X  =  0,1.28,1.92,2.56,  and  2.81  on  February  5. 


Figure  4.30:  Computed  frequency  spectra  of  alongshore  momentum  flux, 
hUV  =  qiU,  at  X  =  0,  1.28,  1.92,  2.56,  and  2.81  on  February  5. 


Figure  4.31:  Computed  frequency  spectra  of  alongshore  bottom  shear  stress  \U\V 
with  fb  =  0.43  at  a;  =  0,  1.28,  1.92,  2.56,  and  2.81  on  February  5. 


Figure  4.32:  Computed  frequency  spectra  of  alongshore  velocity  V  at 
X  =  0, 1.28, 1.92,2.56,  and  2.81  on  February  5  using  linear  (dashed 
curves)  equation  (4.8)  and  nonlinear  (solid  curves)  equation  (3.34). 


Figure  4.34  compares  the  measured  and  computed  cross-shore  variations  of 
the  longshore  current  V  on  February  5  and  6  where  the  overbar  denotes  the  time 
averaging  for  200  <  t  <  500.  The  computed  maximum  and  minimum  alongshore 
velocities,  14iax  and  14un,  during  200  <t<  500  are  also  plotted  in  Figure  4.34  to 
show  that  the  alongshore  velocity  V  is  oscillatory  (positive  and  negative)  outside 
the  surf  zone  and  essentially  unidirectional  (mostly  positive)  in  the  surf  zone  where 
the  temporal  variations  of  V  at  five  different  locations  in  cross-shore  direction  have 
been  shown  in  Figures  4.16  and  4.21.  The  time-dependent  model  with  the  bottom 
friction  factor  /(  =  0.015  calibrated  by  Raubenheimer  et  al.  (1995)  using  runup 
measurements  on  a  natural  beach  also  predicts  the  longshore  current  reasonably 
well.  This  suggests  the  robustness  of  the  adopted  empirical  formula  for  bottom 
friction,  (2.46)  and  (2.47),  although  it  may  be  too  crude  to  predict  the  detailed 
temporal  variation  of  the  bottom  shear  stress  accurately.  Thornton  and  Guza  (1986) 
compared  their  time-averaged  model  with  the  longshore  current  data  and  obtained 
similar  agreement.  No  data  is  available  in  the  swash  zone  to  compare  the  accuracy 
of  the  time-averaged  and  time-dependent  models. 

4.3.5  Computed  Time-Averaged  Quantities 

Figure  4.35  shows  the  computed  cross-shore  variations  of  the  alongshore  ra¬ 
diation  stress  Sxy  given  by  (3.6),  its  cross-shore  gradient,  and  the  time- averaged 
alongshore  bottom  shear  stress  fb\U\V  in  the  time-averaged  alongshore  momen¬ 
tum  equation  (3.5).  The  alongshore  variations  of  the  computed  mean  and  variance 
of  7/  are  found  to  be  negligible  in  (3.5)  relative  to  the  other  two  terms  shown  in 
Figure  4.35.  It  is  noted  that  for  this  2D  model,  the  alongshore  momentum  flux  cor¬ 
rection  n  =  0.  The  assumption  of  alongshore  uniformity  requires  that  the  computed 
mean  and  variance  of  77  at  j/  =  0,  Ay  and  2A7/  should  be  identical.  Figure  4.35  show 
that  the  computed  time-dependent  results  satisfy  the  time-averaged  force  balance 


Figure  4.34:  Comparisons  between  measured  and  computed  longshore  current  V 
together  with  computed  maximum  and  minimum  alongshore  veloci¬ 
ties,  14iax  and  Viain,  on  February  5  and  6. 
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between  the  cross-shore  gradient  of  Sxy  and  the  alongshore  bottom  shear  stress  apart 
from  the  high-frequency  numerical  oscillations  of  dSxyfdx.  The  minute  oscillations 
of  Sxy  caused  by  the  high  frequency  numerical  oscillations  of  the  computed  h,  U, 
and  V  are  almost  indiscernible  in  Figure  4.35  but  are  significantly  amplified  in  the 
cross-shore  derivative  of  Sxy  approximated  by  a  central  finite  difference.  The  along¬ 
shore  bottom  shear  stress  shown  in  Figure  4.35  is  the  maximum  near  the  still  water 
shoreline  as  may  be  expected  from  U  and  V  shown  in  Figures  4.15,  4.20,  4.16  and 
4.21.  Figure  4.35  also  indicates  that  the  numerical  dispersion  of  the  time-dependent 
model  is  negligible  in  the  time-averaged  alongshore  momentum  equation  (3.5)  which 
includes  no  physical  dispersion  or  lateral  mixing  term.  Thornton  and  Guza  (1986) 
already  showed  that  the  lateral  mixing  was  not  important  for  their  longshore  current 
data. 

Figure  4.36  shows  the  computed  cross-shore  variations  of  =  hV  and  h  V. 
The  time-averaged  alongshore  volume  flux  per  unit  width  can  be  expressed  as 

q^  =  hV+{h-h){V  -V)  (4.9) 

where  the  second  term  on  the  right-hand  side  of  (4.9)  is  the  time- averaged  alongshore 
volume  flux  component  due  to  the  oscillatory  components  {h  —  h)  and  [V  —  V). 
Figure  4.36  indicates  that  this  component  in  (4.9)  is  small  in  comparison  to  the 
product  of  the  mean  water  depth  h  and  the  longshore  current  V  except  for  the 
regions  of  small  h  in  the  swash  zone  and  small  V  outside  the  surf  zone. 

To  examine  the  nature  of  large  fluid  velocities  near  the  still  water  shoreline, 
Figures  4.37  and  4.38  show  the  computed  cross-shore  variations  of  the  mean  values 
rf,  U  and  V  as  well  as  the  root-mean-square  values  of  the  oscillatory  components 
(7  ~  7)?  ~  U)  and  {V  —  V)  denoted  by  the  subscript  rms  where  T/rms,  t4ms  and 

14ms  are  the  standard  deviations  of  q,  U  and  V,  respectively.  The  normalized  bot¬ 
tom  elevation  is  indicated  by  the  dotted  line  in  the  top  panel  of  Figures  4.37  and  4.38. 


X 

Figure  4.35:  Computed  cross-shore  variations  of  alongshore  radiation  stress  Sxy, 
its  gradient  and  time-averaged  alongshore  bottom  shear  stress  for 
February  5  and  6. 
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Figure  4.36:  Computed  cross-shore  variations  of  mean  alongshore  volume  flux 
^  =  hV  and  product  of  mean  depth  h  and  longshore  current  V  on 
February  5. 


The  wave  setup  rj  becomes  tangential  to  the  beach  in  the  swash  zone  (Bowen  et  al. 
1968).  The  mean  cross-shore  velocity  U  is  negative  and  represents  the  cross-shore 
return  current  as  explained  by  Kobayashi  et  al.  (1989).  Unlike  the  computed  results 
for  regular  waves  shown  in  Figure  4.5,  T/rms  does  not  decrease  linearly  landward, 
while  Urms  and  I4ns  actually  increase  landward  before  their  decrease  in  the  swash 
zone.  For  lack  of  velocity  data  near  the  still  water  shoreline,  the  computed  cross¬ 
shore  variations  of  t4ms  and  Kms  can  not  be  verified,  whereas  the  cross-shore  one¬ 
dimensional  model  has  been  shown  to  predict  the  free  surface  elevation  t)  in  the 
inner  surf  and  swash  zones  fairly  accurately  (Cox  et  al.  1994;  Raubenheimer  et 
al.  1995).  The  field  data  obtained  at  Torrey  Pines  Beach,  California  by  Guza  and 
Thornton  (1985)  indicated  that  the  variance  was  almost  constant  in  depths 
ranging  from  0.4-2. 5  m  and  that  the  variance  increased  monotonically  as  depth 
decreased  as  shown  in  Figure  4.39,  which  indicated  the  cross-shore  variations  of  the 
variances  of  the  total,  high-frequency,  and  low-frequency  components  of  the  cross¬ 
shore  and  alongshore  velocities.  Their  field  data  and  the  computed  results  shown  in 


Figures  4.37  and  4.38  point  out  the  inadequacy  of  the  conventional  assumption  of 
depth-limited,  linear  breaking  waves  made  in  existing  time-averaged  models. 

The  time-averaged  wave  energy  equation  (3.7)  is  used  to  examine  the  cross¬ 
shore  variations  of  the  time-averaged  values  of  the  specific  energy  the  energy 
flux  the  dissipation  rate  Dj  due  to  bottom  friction,  and  the  dissipation  rate 
Db  due  to  wave  breaking  as  shown  in  Figure  4.40.  The  computed  dissipation  rate 
Db  indicates  intense  wave  breaking  immediately  seaward  of  the  still  water  shoreline 
located  at  x  =  2.56  on  February  5  and  at  x  =  3.82  on  February  6  in  Figure  4.40. 
This  intense  wave  breaking  may  partly  explain  the  large  oscillatory  velocities  near 
the  still  water  shoreline  shown  in  Figures  4.37  and  4.38.  The  computed  wave  energy 
dissipation  in  the  swash  zone  is  caused  mostly  by  bottom  friction  due  to  the  large 
cross-shore  velocity  U  in  the  swash  zone. 

The  large  oscillatory  velocities  near  the  still  water  shoreline  may  also  be 
caused  by  low  frequency  waves.  In  the  very  shallow  water  for  the  velocity  data  of 
Guza  and  Thornton  (1985)  as  shown  in  Figure  4.39,  the  low-  and  high-frequency 
components  of  the  variance  were  of  the  same  magnitude  and  the  variance 
was  dominated  by  the  low  frequency  component.  For  comparison,  the  oscillatory 
components  {'q  —  fj),  {U  —  U)  and  {V  —  V)  are  filtered  to  obtain  the  components 
of  7/^3,  and  in  the  high-frequency  (/'  >  0.05  Hz)  and  low-frequency 

(/'  <  0.05  Hz)  bands  as  shown  in  Figures  4.41  and  4.42.  The  computed  low- 
frequency  components  of  and  in  the  swash  zone  turn  out  to  be  smaller 
than  the  corresponding  high-frequency  components  probably  because  the  incident 
low-frequency  wave  components  are  not  included  in  the  specified  incident  wave  train 
in  the  3  m  depth  as  shown  in  Figure  4.13.  Furthermore,  edge  waves  are  excluded  in 
this  numerical  model.  On  the  other  hand,  the  computed  low-  and  high-frequency 
components  of  1^,3  are  of  the  same  magnitude  near  the  still  water  shoreline  even  in 
the  absence  of  incident  low  frequency  waves  and  edge  waves. 
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Recently,  Kobayashi  et  al.  (1997)  derived  simple  relationships  between  the 
free  surface  and  cross- shore  velocity  statistics  using  linear  progressive  long- wave 
theory.  Accordingly,  the  relationship  between  the  normalized  standard  deviations 
of  77  and  U  may  be  expressed  as 

^  (4- 10) 

whereas  the  mean  cross-shore  velocity  U  is  estimated  as 

(4.11) 

where  T^rms  and  C4ms  are  the  standard  deviation  of  7/  and  U,  respectively,  and  h  is 
the  mean  water  depth  including  wave  setup. 

Figures  4.43  and  4.44  show  the  comparisons  of  the  cross-shore  variations  of 
Unas  and  IJ  computed  using  the  2D  model  and  using  (4.10)  and  (4.11)  together 
with  the  computed  values  of  t/i^s  and  rjrms/h  for  February  5  and  6,  respectively.  The 
normalized  bottom  elevation  is  depicted  with  the  dotted  line  in  the  top  panel  of  these 
figures.  Similar  to  the  experimental  results  of  Kobayashi  et  al.  (1997),  the  computed 
values  of  rjrras/h  in  Figures  4.43  and  4.44  increases  gradually  and  then  rapidly  near 
the  shoreline  as  h  approaches  zero.  Hence,  contrary  to  the  field  data  obtained  by 
Thornton  and  Guza  (1982,  1983),  the  computed  values  of  rj„aslh  do  not  approach 
constant  in  the  inner  surf  zone.  These  figures  also  show  that  the  simple  relationships 
(4.10)  and  (4.11)  based  on  linear  progressive  long-wave  theory  predicts  the  standard 
deviation  of  cross-shore  velocity,  Unas-,  and  the  mean  cross-shore  velocity,  U,  fairly 
well  up  to  the  still  water  shoreline.  Consequently,  the  large  velocities  seaward  of 
the  still  water  shoreline  may  simply  be  explained  by  the  landward  increase  of  rjnas/ h 
in  which  the  landward  decrease  of  rjnas  is  more  gradual  than  that  of  h  as  shown  in 
Figures  4.43  and  4.44.  In  the  swash  zone,  local  nonlinear  effects  such  as  the  bottom 
friction  as  well  as  wave  reflection  appear  to  be  important  in  view  of  the  comparison 
shown  in  these  figures. 
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Figure  4.39:  Total,  high-frequency  and  low-frequency  components  of  cross-shore 
velocity  and  alongshore  velocity  as  a  function  of  mean  depth  h  (from 
Guza  and  Thornton  1985). 


Figure  4.40:  Computed  cross-shore  variations  of  specific  energy  E,  energy  flux 
Ep-,  dissipation  rate  and  Dp  due  to  bottom  friction  and  wave 
breaking,  respectively,  for  February  5  and  February  6. 


Figure  4.41:  Computed  cross-shore  variations  of  variances  of  free  surface  eleva¬ 
tion  ri,  cross-shore  velocity  U  and  alongshore  velocity  V  in  high 
(/'  >  0.05Hz)  and  low  (/'  <  0.05Hz)  frequency  bands  for  Febru- 


Figure  4.42:  Computed  cross-shore  variations  of  variances  of  free  surface  eleva¬ 
tion  7?,  cross-shore  velocity  U  and  alongshore  velocity  V  in  high 
(/'  >  0.05Hz)  and  low  (/'  <  0.05Hz)  frequency  bands  for  Febru- 


Figure  4.43:  Cross-shore  variations  of  /7rms,?/rms/^)  ^rms  and  U  compute^  using 
the  2D  model  as  well  as  cross-shore  variations  of  Unas  and  U  com¬ 
puted  using  (4.10)  and  (4.11)  together  with  the  normalized  bottom 
elevation  for  February  5. 
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Figure  4.44:  Cross-shore  variations  of  7/rms?  and  U  computed  using 

the  2D  model  as  well  as  cross-shore  variations  of  t4ms  and  U  com¬ 
puted  using  (4.10)  and  (4.11)  together  with  the  normalized  bottom 
elevation  for  February  6. 


Chapter  5 


COMPARISON  OF  THE  QUASI-3D  MODEL  WITH 

AVAILABLE  DATA 

The  computed  results  presented  in  Chapter  4  are  based  on  the  depth-integra¬ 
ted  continuity  equation  (3.1)  and  the  depth-integrated  horizontal  momentum  equa¬ 
tions  (3.2)  and  (3.3),  neglecting  the  dispersion  due  to  vertical  non-uniformities  of 
the  horizontal  velocities  u  and  v.  However,  wave  breaking  produces  vertical  vari¬ 
ations  in  the  horizontal  velocities  and  resulting  energy  dissipation  {e.g.,  Svendsen 
and  Madsen  1984).  Furthermore,  the  vertical  variations  of  instantaneous  velocities 
and  shear  stresses  are  required  for  detailed  analyses  of  sediment  transport  in  the  surf 
zone.  Although  Thornton  and  Guza  (1986)  have  already  concluded  that  the  disper¬ 
sion  effects  were  not  important  for  predicting  longshore  currents  on  a  planar  beach 
generated  by  random  waves,  the  computed  results  for  regular  waves  depicted  in 
Figure  4.9  suggest  that  the  dispersion  terms  are  required  to  improve  the  agreement. 

Very  little  is  known  of  these  dispersion  effects  on  surf  zone  hydrodynamics 
apart  from  the  analysis  of  Svendsen  and  Putrevu  (1994)  that  showed  the  importance 
of  the  dispersion  effect  due  to  the  nonlinear  interaction  of  cross-shore  and  longshore 
currents  in  explaining  the  measured  cross-shore  variations  of  longshore  currents 
induced  by  regular  breaking  waves.  They  used  a  time-averaged  model  in  their 
dispersion  analysis  where  linear  wave  theory  with  depth-limited  breaker  height  was 
employed  to  describe  the  wave  motion. 
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To  assess  the  importance  of  the  added  dispersion  terms  m  and  n  in  the  time- 
dependent  momentum  equations  (3.47)  and  (3.48)  due  to  the  vertical  variations  of 
instantaneous  horizontal  velocities  u  and  v,  the  3D  model  described  in  Section  3.2 
is  compared  with  the  same  regular  and  irregular  wave  data  as  in  Chapter  4. 

5.1  Comparison  with  Regular  Wave  Data 

5.1.1  Comparison  of  2D  and  3D  Models  with  Data 

The  3D  model  is  compared  with  the  same  four  experiments  of  Visser  (1991)  as 
in  Section  4.2.  The  depth-averaged  longshore  current  data  of  Visser  (1991)  was  based 
on  the  current  velocity  measurements  of  dye  clouds  injected  at  three  depths.  The 
3D  computations  are  made  in  the  same  way  as  the  corresponding  2D  computations 
presented  in  Section  4.2  except  the  seaward  boundary  location  for  experiment  4 
in  the  2D  and  3D  computations  is  moved  somewhat  closer  to  the  breaker  line  to 
better  satisfy  the  assumption  of  shallow  water.  For  experiment  4,  the  incident  wave 
characteristics  specified  as  input  to  the  numerical  model  are;  water  depth  at  the 
seaward  boundary  d\  =  12.41  cm,  wave  period  T'  =  1.02  s,  wave  height  H'  =  8.4  cm, 
incident  wave  angle  9i  =  13.1°  whereas  the  values  for  experiments  2,  3  and  5  remain 
the  same  as  listed  in  Table  4.1. 

The  temporal  variations  of  the  free  surface  elevation  77  and  the  depth-averaged 
cross-shore  velocity  U  computed  by  the  3D  model  are  practically  the  same  as  those 
computed  by  the  2D  model  shown  in  Figures  4.2  and  4.3  and  hence  are  not  pre¬ 
sented  here.  The  depth-averaged  alongshore  velocity  V  shown  in  Figure  4.4  for  the 
2D  model  is  modified  considerably  as  depicted  in  Figure  5.1  for  experiment  2.  As  a 
whole,  the  alongshore  velocity  V  is  smaller  than  in  Figure  4.4.  The  measured  and 
computed  cross-shore  variations  of  the  local  wave  height  the  wave  setup  57,  and 
the  depth-averaged  longshore  current  V  for  experiments  2-5  are  compared  in  Fig¬ 
ures  5.2  -  5.4,  respectively.  The  normalized  bottom  elevation  for  each  experiment 
is  shown  as  the  straight  line  in  Figure  5.3.  For  the  cross-shore  fiuid  motion  repre- 


sented  by  H  and  rf,  the  differences  between  the  2D  and  3D  computations  are  minor, 
probably  within  the  accuracy  of  these  models  in  light  of  the  agreement  shown  in  Fig¬ 
ures  5.2  and  5.3.  This  implies  that  the  dispersion  term  m  in  the  depth-integrated 
cross-shore  momentum  equation  (3.47)  is  secondary  as  anticipated  by  Kobayashi 
and  Wurjanto  (1992). 

The  average  values  of  the  relative  errors  for  the  2D  and  3D  computations  for 
each  experiment  are  inserted  in  the  parentheses  in  Figures  5.2  and  5.3.  The  relative 
error  for  each  data  point  is  defined  as  \Ym-Yc\l\Ym\  where  Ym  =  measured  value 
and  Yc  =  computed  value.  The  average  relative  error  for  the  mean  water  elevation, 
?7,  is  large  because  of  the  difficulty  in  predicting  the  relatively  small  set-down  and 
setup  accurately. 

On  the  other  hand,  for  the  alongshore  fluid  motion  represented  by  the  long¬ 
shore  current  F,  which  is  dominant  in  comparison  to  the  oscillatory  alongshore 
velocity  in  the  surf  zone  in  these  experiments,  the  dispersion  term  n  in  the  depth- 
integrated  alongshore  momentum  equation  (3.48)  improves  the  agreement  of  the 
longshore  current  profile  noticeably  as  shown  in  Figure  5.4.  The  longshore  currents 
computed  using  fl  =  0.02  and  0.025  for  experiment  4  and  fl  =  0.02  for  experiment  5 
are  also  shown  in  Figure  5.4.  The  alongshore  bottom  shear  stress  Uy  in  (3.48)  is  im¬ 
portant  in  determining  the  magnitude  of  V  but  modifies  its  profile  little  as  expected 
from  the  previous  work  (e.y.,  Longuet-Higgins  1970).  It  is  noted  that  the  cross-shore 
variations  of  H  and  rj  computed  using  these  bottom  friction  factors  for  experiments 
4  and  5  are  indistinguishable  from  those  shown  in  Figures  5.2  and  5.3  for  fl  =  0.05 
and  are  not  plotted  in  these  figures.  The  importance  of  the  dispersion  effect  on 
the  longshore  currents  in  Figure  5.4  is  consistent  with  the  analysis  of  Svendsen  and 
Putrevu  (1994)  based  on  the  vertical  variations  of  currents  only,  while  the  present 
analysis  is  based  on  the  vertical  variations  of  the  instantaneous  horizontal  velocities. 
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Experiment  2 


Figure  5.1:  Computed  temporal  variations  of  depth- averaged  alongshore  velocity 
V  a.t  x  =  0,0.509,0.770,1.550  and  2.265  using  3D  model. 


Data  Computed  (f5=  0.05) 
- 3D: - 2D 


Figure  5.2:  Measured  and  computed  local  wave  height  H  for  four  experiments  for 
comparison  of  2D  and  3D  models. 


118 


Computed  (f^=  0.05) 


X 


Figure  5.4:  Measured  and  computed  depth- averaged  longshore  current  V  for  four 
experiments  for  comparison  of  2D  and  3D  models. 
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To  explain  the  computed  results  in  Figure  5.4,  the  time- averaged  alongshore 
momentum  equation  (3.61)  for  the  case  of  alongshore  uniformity  is  given  by 

(5.1) 

in  which  Sxy  =  hUV  is  the  alongshore  radiation  stress  based  on  the  depth-averaged 
velocities  U  and  V .  The  mean  alongshore  momentum  flux  n  due  to  the  vertical 
variations  of  the  horizontal  velocities  u  and  v  is  included  in  the  3D  model,  whereas 
n  =  0  for  the  2D  model.  Figure  5.5  shows  the  computed  cross-shore  variations  of 
dSxyldx  and  d^jdx  for  experiment  2  as  an  example.  The  cross-shore  gradients  of 
Sxy  for  the  2D  and  3D  computations  are  similar  because  /i,  U  and  V  are  not  affected 
much  by  the  dispersion  terms  m  and  n.  This  figure  indicates  that  the  term  dn/ dx 
included  in  the  3D  model  decreases  the  force  driving  the  longshore  current  near 
the  breaker  point  but  increases  this  force  near  the  shoreline.  Correspondingly,  the 
computed  longshore  current  V  for  the  3D  model  in  Figure  5.4  is  reduced  near  the 
breaker  point  but  increased  near  the  shoreline.  In  short,  the  dispersion  term  n  is 
important  in  the  time-averaged  alongshore  momentum  equation  (5.1). 


5.1.2  Vertical  Variations  of  Longshore  Currents 

The  depth-averaged  longshore  current  data  of  Visser  (1991)  were  based  on 
the  current  velocity  measurements  of  dye  clouds  injected  at  the  surface,  mid-depth, 
and  1  cm  above  the  bed.  The  measured  vertical  variations  of  the  longshore  current  v 
were  presented  for  experiments  3  and  4  in  Figures  5.6  and  5.7  where  the  cross-shore 
location  x  is  indicated  for  each  vertical  measuring  line.  The  vertical  variation  of  the 
alongshore  velocity  v  for  the  3D  model  is  computed  using  (3.58),  (3.59)  and  (3.60)  for 
the  computed  temporal  and  cross-shore  variations  of  F,  and  h  in  these  equations. 
It  is  noted  that  the  2D  model  does  not  yield  the  vertical  velocity  variations.  The 
time- averaged  velocity  u  is  computed  for  the  given  elevation  {z—Z}^)  above  the  bottom 
and  below  the  wave  trough  level  where  zi  =  (— d)  is  the  normalized  bottom  elevation 


introduced  for  clarity.  To  estimate  v  above  the  trough  level,  v  is  also  computed 
for  the  given  value  of  {  =  {z  —  Zb)/h  with  ^  =  1  at  the  instantaneous  free  surface 
and  plotted  as  a  function  of  the  mean  elevation,  {z  —  zi)  =  (h,  with  h  =  mean 
water  depth.  Figures  5.6  and  5.7  show  that  the  vertical  variations  of  v  computed 
for  the  given  [z  —  zi,)  and  are  almost  identical  inside  the  surf  zone  where  the 
breaker  point  was  located  at  a:  =  0.9  for  experiment  3  and  x  =  0.5  for  experiment  4. 
The  measured  and  computed  vertical  profiles  of  v  are  approximately  parallel  except 
outside  the  surf  zone  where  the  assumed  cubic  profile  (3.60)  for  breaking  waves 
may  not  be  appropriate.  Outside  the  surf  zone,  the  vertical  profile  of  v  decreases 
with  the  increase  of  {z  —  Zf,).  In  this  region,  dSxy/dx  =  0,  so  that  the  currents 
are  entirely  driven  by  the  dispersion  mechanism  and  the  (much  weaker)  turbulent 
mixing  (Svendsen  and  Putrevu  1994).  The  discrepancy  between  the  measured  and 
computed  longshore  current  v  inside  the  surf  zone  in  Figures  5.6  and  5.7  is  mostly 
caused  by  the  error  in  the  predicted  depth-averaged  longshore  current  V  shown  in 
Figure  5.4  where  Figure  5.7  is  based  on  fl  =  0.025  for  experiment  4. 
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Figure  5.7:  Measured  and  computed  vertical  variations  for  longshore  current  v  for 
experiment  4  for  3D  model. 


As  for  the  vertical  variation  of  the  instantaneous  cross-shore  velocity  u,  John¬ 
son  et  al.  (1996)  compared  the  numerical  model  based  on  (3.46),  (3.47)  and  (3.52) 
for  normally  incident  waves  with  the  velocity  measurements  below  the  trough  level 
by  Cox  et  al.  (1994).  The  vertical  variations  of  the  measured  and  computed  cross¬ 
shore  velocities  were  relatively  small  below  the  trough  level  and  above  the  bottom 
boundary  layer.  The  computed  vertical  variation  of  u  above  the  trough  level  could 
not  be  verified  for  lack  of  data. 

5.1.3  Numerical  Damping  and  Dissipation 

According  to  Kobayashi  and  Johnson  (1995),  the  value  of  numerical  damping 
coefficient  k  used  for  the  3D  model  in  the  smoothing  procedure  for  breaking  waves 
on  the  gently  sloping  beach  is  on  the  order  of  unity  as  discussed  in  Section  3.2. 
As  a  result,  /c  =  1  has  been  used  unless  stated  otherwise.  The  computation  using 
K  =  1  improves  the  prediction  of  the  transition  zone  but  underpredicts  the  maximum 
setup  for  experiments  4  and  5  in  Figure  5.3.  Figure  5.8  shows  that  the  computed 
wave  setup  rj  using  k  =  0.3  for  the  3D  model  yields  slightly  better  agreement  for 
experiments  4  and  5.  Table  5.1  and  5.2  compares  the  measured  and  computed 
maximum  setup  and  runup  for  the  four  experiments  using  /c=  1  and  0.3  in  the  same 
way  as  in  Table  4.2.  It  is  noted  that  the  measured  values  for  experiment  4  are 
slightly  different  from  those  in  Table  4.2  since  the  normalization  is  made  using  the 
different  value  of  the  wave  height  H'  at  the  shallower  seaward  boundary  location. 
On  the  other  hand.  Figures  5.9  and  5.10  show  that  the  cross-shore  variations  of  the 
wave  height  H  and  the  depth-averaged  longshore  current  V  based  on  k  =  0.3  do 
not  change  much  as  compared  to  the  computed  results  with  k  =  1.  In  summary, 
the  numerical  damping  coefficient  k,  in  the  range  of  k  =0.3-1  modifies  the  computed 
results  slightly  but  does  not  change  the  essential  features  of  the  computed  results 
using  the  3D  model. 


Table  5.1:  Measured  and  computed  maximum  setup  rf  using  different  values  of  n 
for  3D  model. 
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Table  5.2:  Measured  and  computed  runup  using  different  values  of  /c  for  3D  model. 
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The  time-averaged  energy  quantities  are  shown  in  Figure  5.11.  The  numerical 
dissipation  rate,  Dnumericai,  shown  in  the  last  panel  and  defined  in  (3.69)  is  the  differ¬ 
ence  between  the  values  of  Db  computed  using  the  time-averaged  energy  equation 
(3.63)  and  the  physical  dissipation  rate  Db  predicted  using  (3.68)  for  the  3D  model. 
It  is  noted  that  the  numerical  dissipation  rate  is  of  the  same  order  as  the  wave 
breaking  dissipation  rate  computed  using  (3.68)  explicitly  in  the  3D  model.  The 
values  of  K  =  0.3  and  1.0  have  only  minor  effects  on  the  time- averaged  energy  quan¬ 
tities.  This  suggests  the  shortcoming  of  the  cubic  velocity  profile  assumed  in  (3.60) 
to  describe  the  wave  energy  dissipation  due  to  breaking  on  the  gentle  slope.  This 
probably  arises  from  the  fact  that  the  wave  front  (roller)  is  not  modeled  specifically 
in  the  3D  model  (Johnson  et  al.  1996). 


Figure  5.8:  Measured  and  computed  wave  setup  t)  for  four  experiments  using  dif 
ferent  values  of  k  for  3D  model. 


Figure  5.11:  Computed  cross-shore  variations  of  time- averaged  specific  energy  E, 
energy  flux  Ep^  bottom  frictional  dissipation  rate  D/,  breaking  dis¬ 
sipation  rate  Db  predicted  physically  by  3D  model,  and  numeri¬ 
cal  dissipation  rate  ^numerical  estimated  using  time-averaged  energy 
equation. 


5.1.4  Computed  Instantaneous  Velocity  Field 

Figure  5.12  shows  the  computed  cross-shore  variations  of  the  free  surface 
77,  the  depth- averaged  cross-shore  velocity  U,  the  near-bottom  cross-shore  velocity 
correction  Ub,  the  depth- averaged  alongshore  velocity  V,  and  the  near-bottom  along¬ 
shore  velocity  correction  Vb  at  five  time  levels  throughout  the  final  wave  period  for 
experiment  2.  The  computed  cross-shore  variations  at  t  =  299  and  t  =  300  are 
identical  due  to  periodicity.  The  saw  tooth  profile  develops  as  the  wave  propagates 
shoreward.  The  near-bottom  cross-shore  velocity  correction,  Ub,  is  assumed  to  be 
out  of  phase  with  U  in  (3.54)  to  ensure  that  |ui,|  <  \U\  where  Ub  =  {U  +  Ub)  is 
the  near  bottom  cross-shore  velocity.  Consequently,  the  value  of  Ub  computed  using 
(3.54)  changes  abruptly  as  U  changes  its  sign.  The  magnitude  of  Vb  is  small  in 
comparison  to  the  magnitude  of  V,  suggesting  that  the  deviation  of  the  alongshore 
velocity  v  given  by  (3.58)  from  the  depth-averaged  velocity  V  is  fairly  small  as  is 
the  case  with  the  longshore  current  v  shown  in  Figures  5.6  and  5.7. 

5.2  Comparisons  with  Irregular  Wave  Data 

The  3D  computations  for  the  Leadbetter  beach  data  of  Thornton  and  Guza 
on  February  5  and  6  are  made  in  the  same  way  as  the  corresponding  2D  computations 
in  Section  4.3.  The  bottom  friction  factor  is  taken  to  be  fl  =  0.015  for  the  2D  and 
3D  computations.  The  3D  computations  are  also  made  using  fl  =  0.01. 

Figure  5.13  compares  the  measured  and  computed  cross-shore  variations  of 
the  normalized  root-mean-square  wave  height  Hrms  on  February  5  and  6.  The  differ¬ 
ences  between  the  2D  and  3D  computations  are  less  than  those  shown  in  Figure  5.2 
for  the  regular  waves.  The  cross-shore  variations  of  iTrms  computed  using  fl  =  0.01 
are  indistinguishable  from  those  shown  in  Figure  5.13  for  fl  =  0.015. 

Figure  5.14  shows  the  computed  depth-averaged  longshore  currents  V  using 
fl  =  0.01  and  0.015  in  comparison  with  the  longshore  current  measured  at  a  distance 
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Figure  5.12:  Computed  cross-shore  variations  of  7,  iii,  V,  and  Vh  at  5  time  levels, 
t  =  299.0,  299.25,  299.5,  299.75,  and  300.0  for  experiment  2. 


above  the  bottom.  These  comparisons  are  meaningful  only  if  the  longshore  current 
V  varies  little  vertically.  Figure  5.15  shows  the  computed  vertical  variations  of  v 
based  on  fl  =  0.01  at  six  cross-shore  locations  on  February  5  where  x  =  2.56  and 
2.81  are  in  the  swash  zone.  It  is  noted  that  the  time-averaged  value  of  v  at  the  given 
distance  [z  —  Zh)  from  the  bottom  is  computed  only  if  this  elevation  is  wet  always 
during  200  <  t  <  500,  whereas  the  time-averaged  value  v  at  the  given  normalized 
distance  (  =  [z  —  Zb)fh  is  computed  as  long  as  h  >  0  and  plotted  at  the  mean 
elevation,  (z  —  Zb)  =  (h,  where  u  =  0  set  during  h  =  0  is  included  in  the  time 
averaging.  The  computed  vertical  variations  of  v  on  February  6  are  also  small  and 
the  measured  longshore  currents  may  be  assumed  to  represent  the  depth-averaged 
longshore  currents.  It  is  noted  that  the  vertical  variations  of  v  for  irregular  waves 
in  Figure  5.15  are  smaller  than  those  for  regular  waves  in  Figures  5.6  and  5.7. 

Figure  5.14  shows  that  the  dispersion  effects  included  in  the  3D  model  im¬ 
prove  the  agreement  somewhat  if  the  bottom  friction  factor  fl  is  reduced  to  fl  =  0.01 
from  fl  =  0.015  used  for  the  2D  model.  Moreover,  the  dispersion  effects  on  the  long¬ 
shore  currents  induced  by  breaking  irregular  waves  are  secondary  in  comparison  to 
those  shown  in  Figure  5.4  for  regular  waves.  To  confirm  this  conclusion.  Figure  5.16 
shows  the  computed  cross-shore  variations  of  dS^y/dx  and  dn/dx  on  February  5 
in  the  same  way  as  Figure  5.5  for  regular  waves.  The  term  dnjdx  included  in  the 
3D  model  is  secondary  in  the  alongshore  momentum  equation  (5.1).  Thornton  and 
Guza  (1986)  already  showed  that  the  dispersion  or  lateral  mixing  was  not  impor¬ 
tant  for  their  longshore  current  data  where  the  lateral  mixing  term  was  based  on 
the  eddy  viscosity  formulation  of  Longuet-Higgins  (1970). 

Figure  5.17  shows  the  time-averaged  energy  quantities  for  February  5.  Similar 
to  Figure  5.11  for  regular  waves,  the  numerical  dissipation  rate  in  the  3D  model  is 
appreciable  relative  to  the  physical  dissipation  rates  Dj  and  Db  for  irregular  waves 


Figure  5.14:  Measured  and  computed  depth- averaged  longshore  current  V  on 
February  5  and  6  for  comparison  of  2D  and  3D  models. 


X 


Figure  5.17:  Computed  cross-shore  variations  of  time-averaged  specific  energy  E, 
energy  flux  Ep,  bottom  frictional  dissipation  rate  Dj,  wave  breaking 
dissipation  rate  Dp  predicted  physically  by  3D  model,  and  numeri¬ 
cal  dissipation  rate  ^numerical  estimated  using  time-averaged  energy 
equation  for  February  5. 
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Chapter  6 


COMPARISON  OF  THE  3D  MODEL  WITH  BARRED 

BEACH  DATA 


The  concept  of  alongshore  radiation  stress  has  been  shown  to  satisfactorily 
describe  the  cross-shore  distribution  of  longshore  current  on  a  planar  beach  (e.^., 
Svendsen  and  Putrevu  1994).  The  peak  is  located  near  the  breaking  zone  inside 
the  surf  zone  where  the  cross-shore  gradient  of  wave  height  is  maximum.  For  regu¬ 
lar  waves,  horizontal  mixing  is  required  for  smoothing  the  current  profile,  whereas 
Thornton  and  Guza  (1986)  found  that  for  irregular  waves,  horizontal  mixing  is  not 
important  as  a  result  of  the  randomness  of  the  breaker  locations. 

On  a  barred  beach,  conceptually,  waves  will  break  on  the  bar,  reform  and 
break  again  on  the  beach  face  producing  two  peaks  in  the  longshore  current  distri¬ 
bution.  Contrary  to  this  concept,  the  measurements  of  longshore  currents  on  the 
barred  beach  obtained  during  the  Duck  Experiment  on  Low-Frequency  and  Incident- 
Band  Longshore  and  Across-shore  Hydrodynamics  (delilah)  experiment  (Smith  et 
al.  1993)  generally  indicated  a  broad  peak  in  the  bar  trough  region.  Existing  time- 
averaged  models  for  longshore  currents,  which  couple  four  governing  equations  for 
the  wave  height,  wave  angle,  mean  water  surface  elevation,  and  longshore  current, 
have  not  been  able  to  predict  these  broad  peak  longshore  current  data  (Smith  et  al. 
1993). 

Smith  et  al.  (1993)  developed  a  one- dimensional  time-averaged  numerical 
model  for  longshore  current  that  included  the  effect  of  turbulence  due  to  wave 


breaking  through  a  general  transport  equation  for  the  nnean  turbulent  kinetic  en¬ 
ergy.  Their  model  improved  the  model  proposed  by  Larson  and  Kraus  (1991)  but 
produced  an  unrealistic  high  peak  on  the  beach  face.  Church  and  Thornton  (1993) 
developed  a  model  using  a  spatially  varying  bottom  friction  coefficient  based  on 
a  one-dimensional  turbulent  kinetic  energy  equation  associated  with  the  breaking- 
wave  induced  turbulence.  However,  this  model  was  unable  to  satisfactorily  predict 
the  broad  peak  of  the  longshore  current  distribution  observed  in  the  DELILAH  exper¬ 
iment.  Momentum  fluxes  associated  with  mass  transport  above  the  trough  level  of 
broken  waves,  which  were  ignored  in  the  other  models,  were  included  in  the  model 
developed  by  Kuriyama  (1994).  His  model  with  additional  empirical  coefficients 
was  compared  with  the  longshore  currents  measured  at  Hazaki  Oceanographical 
Research  Facility  (hoRF)  in  Japan.  At  present,  there  is  no  model  available  to  pre¬ 
dict  the  broad  peak  of  the  longshore  current  on  a  barred  beach  in  a  physically 
satisfactory  manner. 

To  assess  whether  the  time-dependent  3D  numerical  model  including  the 
dispersion  effects  in  Section  3.2  is  capable  of  predicting  the  longshore  current  on  a 
barred  beach,  the  3D  model  is  compared  with  the  DELILAH  field  data  of  Smith  et 
al.  (1993).  The  comparison  suggests  that  the  broad  peak  of  the  longshore  current 
in  the  bar  trough  region  cannot  be  explained  under  the  assumption  of  alongshore 
uniformity. 

6.1  Estimation  of  Incident  Irregular  Waves 

The  3D  model  is  compared  with  the  DELILAH  data  on  October  14,  1990 
at  1900  EST,  which  included  the  cross-shore  variations  of  the  measured  root-mean- 
square  wave  height  and  longshore  current.  The  frequency  spectrum  measured  at  the 
8  m  water  depth  was  narrow  banded  in  frequency  with  symmetric  directional  distri¬ 
butions  about  a  mean  oblique  wave  direction.  The  wave  conditions  at  the  8  m  depth 


were:  the  root-mean-square  wave  height  =  0.83  m;  the  spectral  peak  period 
Tp  =  12.0  sec;  and  the  dominant  incident  wave  direction  Oi  =  18°.  The  bathymetry 
was  nearly  uniform  in  the  alongshore  direction. 

The  seaward  boundary  of  the  numerical  model  based  on  the  assumption  of 
shallow  water  waves  is  taken  at  the  water  depth  d'  =  3.64  m  below  the  still  water 
level  where  the  measured  root-mean-square  wave  height  was  1.02  m.  The  mea¬ 
sured  frequency  spectrum  at  d'  =  8  m  is  used  to  estimate  the  assumed  unidirectional 
frequency  spectrum  and  the  predominant  incident  wave  direction  at  d'  =  3.64  m  us¬ 
ing  the  computer  program  RESHOAL  developed  by  Poff  and  Kobayashi  (1993)  as 
explained  in  the  following. 

RESHOAL  assumes  a  straight  shoreline  with  parallel  bottom  contours.  For  a 
given  incident  directional  random  wave  spectrum  at  a  deeper  water  depth,  RESHOAL 
computes  the  directional  random  wave  spectrum  at  a  specified  shallow  water  depth 
using  linear  finite-depth  wave  theory  for  directional  random  wave  shoaling  and  re¬ 
fraction  (LeMehaute  and  Wang  1982).  The  incident  directional  random  wave  spec¬ 
trum  at  the  deeper  water  depth  d'  =  8  m  is  assumed  to  be  given  by  the  product  of 
the  TMA  frequency  spectrum  (Bouws  et  al.  1985)  and  the  Mitsuyasu-type  direc¬ 
tional  spreading  function  (Goda  1985).  The  input  parameters  for  RESHOAL  at  the 
deeper  water  depth  d'  =  8  m  are:  =  spectral  estimate  of  significant  wave  height; 

Tp  —  spectral  peak  period  {T^  =  12  s  for  this  data);  7  =  spectral  peak  enhance¬ 
ment  factor;  6i  =  dominant  incident  wave  direction  {9i  =  18°);  Smax  =  maximum 
value  of  the  spreading  parameter.  RESHOAL  computes  the  directional  spectrum, 
frequency  spectrum  and  directional  spreading  function  at  the  shallower  water  depth 
d'  =  3.64  m.  The  parameters  7,  and  s^ax  need  to  be  calibrated  such  that 
the  root-mean-square  of  wave  height  at  the  3.64  m  depth  is  equal  to  the  measured 
value  of  =1-02  m  and  the  assumed  incident  directional  wave  spectrum  at  the 
8  m  water  depth  is  similar  to  the  measured  spectrum.  The  calibrated  values  are 


=  1.35  m;  7  =5;  and  5max=120  at  the  8  m  water  as  shown  in  Figure  6.1. 
It  is  noted  that  the  assumption  of  =  H^^(y/2  yields  —  0.95  m  for 
=  1.35  m,  which  is  slightly  larger  than  the  measured  value  =  0.83  m. 

Figures  6.1  shows  the  measured  and  fitted  frequency  and  directional  spectra 
at  the  8  m  water  depth  and  the  computed  frequency  and  directional  spectra  at 
the  seaward  boundary  d'  =  3.64  m.  The  fitted  directional  spectrum  at  the  8  m 
depth  is  the  TMA  frequency  spectrum  with  7  =  5  and  =  1.35  m  with  the 
Mitsuyasu-type  directional  spreading  function  with  Smax  =  120.  The  computed 
frequency  spectrum  with  =  11.9  s  at  d'  =  3.64  m  is  used  to  compute  the  incident 
wave  trains  at  the  seaward  boundary  using  (4.3)  required  as  the  input  to  the  3D 
model.  This  incident  frequency  wave  spectrum  does  not  include  low-frequency  wave 
components  as  shown  in  Figure  6.1.  The  computed  dominant  incident  wave  direction 
is  6i  —  12°  at  d'  =  3.64  m  as  may  be  seen  from  the  computed  directional  spectrum 
at  d'.  =  3.64  m  in  Figure  6.1  which  suggests  that  the  assumption  of  unidirectional 
random  waves  may  be  reasonable. 

Similar  to  the  computations  made  in  Chapters  4  and  5,  the  normalized  com¬ 
putation  duration  is  taken  as  tmax  =  500  corresponding  to  t^^ax  —  99  min.  The 
sampling  rate  At'  in  (4.3)  is  taken  to  be  the  same  as  the  sampling  rate  of  the 
field  data,  At'  =  0.125  s.  The  bottom  friction  factor  is  assumed  to  be  /j  =0.015. 
The  computed  results  presented  in  the  following  section  are  based  on  the  normal¬ 
ization  using  the  wave  conditions  at  the  seaward  boundary  of  d'  =  3.64  m,  i.t., 
the  measured  root-mean-square  wave  height  H'  =  =  1.02  m;  the  computed 

spectral  peak  period  T'  =  Tp=11.9  sec;  and  the  computed  dominant  wave  direc¬ 
tion  9i  =  12°.  Correspondingly,  cr  =  T'{glH')2  =  37  and  Be  =  Bi  —  0.21  in 
radians.  The  assumptions  of  >>  1  and  <C  1  are  satisfied  for  this  data. 
The  normalized  grid  spacings  are  taken  as  Ax  ~  Ay  =  0.0106  corresponding  to 
the  dimensional  cross-shore  and  alongshore  grid  spacings  of  Ax'  =  0.40  m  and 
Ay'  =  1.91  m,  respectively. 


Figure  6.1;  Measured  and  fitted  frequency  and  directional  spectra  at  the  deeper 
depth  d'  —  Sm  and  the  shoaled  and  refracted  frequency  and  directional 
spectra  at  the  shallower  depth  d'  =  3.64  m. 
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6.2  Comparison  with  Measurements 

Figure  6.2  shows  the  comparisons  between  the  measured  and  computed  cross¬ 
shore  variations  of  the  local  root-mean-square  wave  height  i/rms  and  the  longshore 
current  V  together  with  the  measured  bottom  profile.  The  computed  temporal 
variation  of  77  for  200  <  t  <  500  is  used  to  obtain  Urms  based  on  the  zero-up  crossing 
method  whereas  the  computed  longshore  current  V  is  obtained  by  averaging  the 
temporal  variation  of  the  depth-averaged  alongshore  velocity  V  for  the  duration 
200  <  t  <  500. 

Figure  6.2  shows  that  the  3D  model  without  the  incident  low-frequency  wave 
components  underpredicts  the  root-mean-square  wave  height  in  the  bar  trough  re¬ 
gion.  Moreover,  the  model  predicts  a  peak  in  the  longshore  current  at  the  seaward 
edge  of  the  bar  crest  in  contrast  to  a  broad  peak  in  the  bar  trough  region. 

The  computed  cross-shore  variations  of  the  time-averaged  energy  quantities 
in  (3.63)  are  shown  in  Figure  6.3.  The  dissipation  due  to  bottom  friction  £>/  has 
a  peak  near  the  shoreline  where  the  oscillatory  cross-shore  velocity  is  large.  The 
dissipation  rate  due  to  wave  breaking,  Bjb,  shown  in  the  fourth  panel  has  peaks 
near  the  bar  crest  and  on  the  beach  face  where  the  waves  break.  The  numerical 
dissipation  rate  shown  in  the  last  panel  is  dominant  over  the  physical  dissipation 
rates  computed  explicitly  in  the  3D  model  and  is  the  maximum  at  the  seaward  edge 
of  the  bar  crest.  This  numerical  dissipation  rate  appears  to  compensate  for  the 
underpredicted  value  of  Ds  at  the  bar  crest. 

Figure  6.4  shows  the  computed  cross-shore  variations  of  dSxyldx  and  dnfdx  in 
the  alongshore  momentum  equation  (5.1).  The  dispersion  term  dnjdx  is  secondary 
in  comparison  to  the  cross-shore  gradient  of  the  alongshore  radiation  stress,  dS^yfdx. 
The  driving  force  represented  by  dSxyj dx  is  small  in  the  bar  trough  region. 

As  an  attempt  to  explain  the  broad  peak  in  the  longshore  current  distribution, 
the  effects  of  incident  low-frequency  waves  and  alongshore  non-uniformity  of  incident 


wind  waves  on  the  cross-shore  distribution  of  longshore  current  are  examined  in  the 
subsequent  sections. 

6.3  Effects  of  Incident  Low-Frequency  Waves 

The  incident  wave  spectrum  at  the  seaward  boundary  shown  in  Figure  6.1 
does  not  include  low-frequency  components.  The  low-frequency  components  might 
modify  wave  breaking  on  the  bar  crest  and  improve  the  agreement  of  the  longshore 
current  profile.  As  a  first  approximation,  uniform  low-frequency  components  are 
added  to  the  incident  wave  spectrum  as  shown  in  Figure  6.5.  The  additional  energy 
of  these  low-frequency  components  is  taken  to  be  4%  and  20%  of  the  energy  of 
the  wind  wave  frequency  components  to  examine  the  sensitivity  of  the  computed 
results  to  the  amount  of  the  low-frequency  components.  It  is  noted  that  the  20% 
low-frequency  component  at  the  3.64  m  depth  are  expected  to  be  too  large  in  reality. 

The  measured  and  computed  cross-shore  variations  of  the  root-mean-square 
wave  height  and  longshore  current  are  compared  in  Figure  6.6.  The  agreement  for 
the  root-mean-square  wave  height  is  improved  somewhat  but  the  additional  low- 
frequency  components  modify  the  longshore  current  profile  little.  Consequently, 
the  broad  peak  of  the  longshore  current  in  the  bar  trough  cannot  be  explained  by 
incident  low-frequency  waves. 

6.4  Effects  of  Alongshore  Non-Uniformity  of  Incident  Wind  Waves 

Longshore  currents  have  been  primarily  modeled  assuming  alongshore  unifor¬ 
mity,  although  it  has  been  known  that  alongshore  non-uniformities  affect  longshore 
currents.  Putrevu  et  al.  (1995)  derived  a  semi-analytical  solution  for  longshore 
current  that  allows  for  weak  alongshore  variation  in  the  bottom  topography.  Their 
calculation  showed  that  the  longshore  current  could  deviate  by  up  to  30%  from  the 
mean  for  a  10%  deviation  of  the  bottom  topography.  Symonds  and  Huntley  (1980) 
showed  that  the  alongshore  gradient  of  wave  setup  can  shift  the  peak  of  the  long- 
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Figure  6.2:  Measured  and  computed  cross-shore  variations  of  root-mean-square 
wave  height  J^rms  and  longshore  current  V  together  with  normalized 
bottom  profile  for  DELILAH  experiment. 
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Figure  6.3:  Computed  cross-shore  variations  of  time- averaged  specific  energy  E, 
energy  flux  Ep,  bottom  frictional  dissipation  rate  £>/,  wave  breaking 
dissipation  rate  Dp  predicted  physically  by  3D  model,  and  numer¬ 
ical  dissipation  rate  D„umericai  estimated  using  time-averaged  energy 
equation. 
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Figure  6.5:  Incident  wave  frequency  spectrum  at  the  seaward  boundary  including 
the  4%  and  20%  low-frequency  components. 
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shore  current  shoreward.  Gourlay  (1976),  Keeley  and  Bowen  (1977),  and  Wu  et 
al.  (1985)  studied  the  effect  of  the  alongshore  variation  of  breaker  height  on  the 
longshore  current.  Gourlay  (1976)  conducted  experiments  in  a  laboratory  to  study 
a  non-uniform  longshore  current  system  generated  by  the  alongshore  gradient  of 
breaker  height  behind  an  offshore  breakwater.  Keeley  and  Bowen  (1977)  measured 
the  longshore  currents  in  Martinique  Beach,  Canada  at  several  alongshore  locations 
and  found  that  alongshore  variation  in  the  breaker  height  contributed  about  10% 
to  the  overall  current  strength.  Wu  et  al.  (1985)  compared  their  two-dimensional 
circulation  model  with  the  Nearshore  Sediment  Transport  Study  (NSTS)  field  mea¬ 
surements  at  Leadbetter  Beach,  California  and  demonstrated  that  the  non-uniform 
wave  field  resulted  in  the  formation  of  a  non-uniform  longshore  current  pattern. 

To  study  the  effect  of  the  alongshore  variation  of  incident  wave  conditions 
within  the  limitation  of  the  3D  model  based  on  three  cross-shore  lines  shown  in 
Figure  3.1,  the  incident  wave  trains  rn  specified  as  new  input  to  the  model  are 
modified  as  follows: 

Line  1:  new  (77^)1  =  old  (771)1 

Line  2:  new  (771)2  =  old  (771)2  X  (1  —  S^,) 

Line  3:  new  (771)3  =  old  (771)3  x  (1  —  26^,) 

where  the  old  time  series  (77^)1,  (771)2  and  (77^)3  have  been  computed  using  (4.3)  for  in¬ 
cident  unidirectional  random  waves  of  alongshore  uniformity.  The  distance  between 
two  adjacent  lines  is  Ay'  =  1.91  m.  The  dimensionless  parameter  Sr,  is  taken  to  be 
much  less  than  unity  to  satisfy  the  assumption  of  gradual  alongshore  variation.  The 
incident  wave  intensity  decreases  or  increases  in  the  down- wave  direction  depending 
on  >  0  or  6,,  <  0,  respectively. 

The  computed  results  using  these  new  incident  wave  trains  are  shown  in 

Figures  6.7  and  6.8  for  the  cases  of  6,,  =  0,  0.0005  and  0.001  and  for  the  cases  of 

(5^  =  0  and  -0.001,  respectively.  The  root-mean-square  wave  height  changes  very 
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little  since  the  specified  change  in  the  incident  wave  train  is  very  small.  For  the 
wave  intensity  decreasing  in  the  down- wave  direction,  the  longshore  current  profile 
increases  almost  uniformly  across  the  shoaling  region  and  over  the  bar  crest.  The 
increase  in  the  longshore  current  is  larger  in  the  bar  trough  region,  whereas  the 
increase  is  smaller  in  the  swash  zone.  The  broad  peak  in  the  bar  trough  region  is 
similar  to  the  broad  peak  observed  in  the  field.  For  the  wave  intensity  increasing 
in  the  down-wave  direction  as  shown  in  Figure  6.8  the  longshore  current  in  the  bar 


trough  region  is  decreased  significantly  and  becomes  negative. 

To  explain  the  computed  results  in  Figures  6.7  and  6.8,  the  time-averaged 
alongshore  momentum  equation  (3.61)  is  rewritten  as 


—S 


dn  -^dTj 


1  d  - -r2  _ 


in  which  Sxy  —  hUV  is  the  alongshore  radiation  stress  based  on  the  depth-averaged 
velocities  U  and  V.  Figure  6.9  shows  the  cross- shore  variations  of  the  driving  forces 
on  the  left  hand  side  of  (6.1)  for  the  cases  of  6r,  =  0.001  and  —0.001.  The  terms  of 
dn/dx  and  |5(r7  -  fj^ldy  are  on  the  order  of  0.005  or  less  in  the  bar  trough  region 
and  secondary  in  comparison  to  the  other  two  terms  plotted  in  this  figure.  For  both 
cases,  the  cross-shore  gradient  of  the  alongshore  radiation  stress,  dSxy/dx,  driving 
the  longshore  current  is  very  small  in  the  bar  trough  region.  The  additional  term 
hdrjfdy  associated  with  the  alongshore  wave  setup  gradient  modifies  the  driving 
force  significantly  in  the  bar  trough  region.  The  incident  wave  intensity  and  resulting 
wave  setup  decrease  or  increase  in  the  down  wave  direction  depending  on  >  0  or 


(5.,,  <  0,  respectively. 

Figures  6.10  and  6.11  show  the  computed  vertical  variations  of  the  long¬ 
shore  current  v  for  the  alongshore  uniform  case  =  0)  and  the  non-uniform  case 
(,5^,j  =  0.001),  respectively,  computed  in  the  same  way  as  for  the  planar  beach 

case  discussed  in  Section  5.2.  Similar  to  the  computed  results  for  irregular  waves 
on  the  planar  beach  shown  in  Figure  5.15,  the  computed  vertical  variations  of  v  are 


very  small.  However,  the  longshore  currents  measured  at  the  same  barred  beach 
by  Haines  and  Sallenger  (1994)  showed  marked  vertical  variations,  which  cannot  be 
explained  by  the  present  3D  model. 

To  examine  the  effects  of  the  alongshore  non-uniformity  on  planar  beaches, 
the  modified  incident  wave  trains  for  the  cases  of  =  0.0005  and  0.001  are  also 
specified  for  the  computations  for  the  regular  wave  experiment  2  of  Visser  (1991) 
and  the  irregular  wave  data  of  Thornton  and  Guza  (1986)  on  February  5.  For  the 
planar  beach  as  shown  in  Figure  6.12,  the  longshore  current  increases  almost  uni¬ 
formly  in  the  shoaling  and  surf  zones  except  in  the  swash  zone.  Contrary  to  the 
computed  results  for  the  barred  beach  shown  in  Figure  6.7,  the  longshore  current 
profile  shape  on  the  planar  beaches  is  not  sensitive  to  the  alongshore  non-uniformity. 
To  explain  this  difference.  Figure  6.13  shows  the  cross-shore  variations  of  the  driving 
forces  in  the  time-averaged  momentum  equation  (6.1)  for  both  regular  and  irregu¬ 
lar  waves  on  the  planar  beaches.  The  additional  driving  force  terms  hdfjidy  and 
_  rfyidy  due  to  the  alongshore  variations  are  small  in  comparison  to  the  main 
driving  force,  dS^yldx,  causing  the  almost  uniform  increase  in  the  longshore  current 
without  changing  its  shape  as  shown  in  Figure  6.12.  The  computed  results  discussed 
above  imply  that  the  broad  peak  of  the  longshore  current  on  a  barred  beach  can  be 
caused  by  the  very  small  alongshore  variation  of  wave  height  and  setup.  This  may 
explain  why  existing  longshore  current  models  based  on  the  assumption  of  along¬ 
shore  uniformity  were  regarded  to  be  adequate  before  their  comparisons  with  the 
barred  beach  data.  For  planar  beaches,  the  effect  of  alongshore  non-uniformity,  even 
if  it  exists,  can  be  accounted  for  by  adjusting  the  constant  bottom  friction  factor 
that  does  not  change  the  longshore  current  profile  shape  as  shown  in  Figures  5.4  and 
5.14.  On  the  other  hand,  for  barred  beaches,  the  very  small  alongshore  variation  of 
wave  height  and  setup  modifies  the  longshore  current  profile  shape  that  cannot  be 
changed  much  by  adjusting  the  constant  bottom  friction  factor. 


Figure  6.6:  Effects  of  low-frequency  wave  components  on  the  cross-shore  variations 
of  the  root-mean-square  wave  height  Hyms  and  longshore  current  V  for 
barred  beach. 


Figure  6.7:  Effects  of  alongshore  decrease  of  incident  wave  intensity  on  cross-shore 
variations  of  root-mean-square  wave  height  Hims  and  longshore  current 
V  for  barred  beach. 


Figure  6.8:  Effects  of  alongshore  increase  of  incident  wave  intensity  on  cross-shore 
variations  of  root-mean-square  wave  height  ifrms  and  longshore  current 
V  for  barred  beach. 
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Figure  6.9:  Computed  cross-shore  variations  of  driving  forces  in  time-averaged 
alongshore  momentum  equation  (6.1)  for  6,,  =  0.001  with  drj/dy  <  0 
(wave  setup  decreasing  in  down-wave  direction)  and  for  ^.^,=-0.001  with 
dfildy  >  0. 


Computed  vertical  variations  of  longshore  current  v  for  alongshore 
uniform  case  =  0)  for  barred  beach. 
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Figure  6.11:  Computed  vertical  variations  of  longshore  current  v  for  alongshore 
non-uniform  case  =  0.001)  for  barred  beach. 


Figure  6.12:  Effects  of  alongshore  decrease  of  incident  wave  intensity  on  the  Ion; 
shore  current  profile  for  planar  beaches. 


Chapter  7 


SUMMARY  AND  CONCLUSIONS 


The  finite-amplitude  shallow-water  equations  under  the  assumption  of  small 
incident  angles  are  solved  numerically  to  predict  the  cross-shore  and  temporal  vari¬ 
ations  of  the  free  surface  elevation  and  the  cross-shore  and  alongshore  velocities  in 
the  swaish  and  surf  zones.  The  use  of  the  finite-amplitude  shallow-water  equations 
limits  the  computation  domain  within  a  short  distance  from  the  shoreline.  The 
assumption  of  small  incident  angles  reduces  computational  efforts  considerably  and 
avoids  difl&culties  associated  with  lateral  boundary  conditions.  This  assumption 
may  be  too  restrictive  for  field  applications  but  allows  us  to  compute  the  dominant 
cross-shore  fluid  motion  along  each  cross-shore  line  using  the  one-dimensional  con¬ 
tinuity  and  cross-shore  momentum  equations  for  normally  incident  waves  developed 
by  Kobayashi  et  al.  (1987,  1989). 

The  secondary  alongshore  velocity,  which  may  vary  slowly  in  the  alongshore 
direction,  is  then  computed  using  the  alongshore  momentum  equation  for  the  com¬ 
puted  free  surface  elevation  and  cross-shore  velocity  along  the  cross-shore  lines.  The 
time-averaged  alongshore  momentum  equation  is  used  to  check  the  accuracy  of  the 
numerical  method  as  well  as  the  seaward  and  landward  boundary  algorithms  em¬ 
ployed  to  solve  the  time-dependent  alongshore  momentum  equation. 

Two  models  are  developed  in  this  study.  The  2D  model  neglects  the  vertical 
variations  of  the  instantaneous  horizontal  velocities,  whereas  the  quasi-3D  model 
includes  the  dispersion  effects  due  to  the  vertical  variations  of  the  horizontal  veloci- 


ties  on  the  depth-integrated  cross-shore  and  alongshore  momentum  equations  (3.47) 
and  (3.48).  The  dispersion  terms  m  and  n  in  these  equations  express  the  additional 
cross-shore  and  alongshore  momentum  fluxes,  respectively.  It  is  noted  that  for  the 
2D  model,  m  =  0  and  n  =  0.  Two  additional  equations  for  m  and  n  are  derived  from 
the  depth-dependent  cross-shore  and  alongshore  momentum  equations.  To  close  the 
problem,  the  deviations  of  the  horizontal  velocities  from  the  depth-averaged  veloci¬ 
ties  are  assumed  to  be  expressed  by  the  cubic  profile  normalized  by  the  near-bottom 
deviations  and  the  instantaneous  water  depth.  The  equations  for  m  and  n  are  solved 
numerically  together  with  the  cross-shore  and  alongshore  momentum  equations  to 
compute  the  temporal  and  cross-shore  variations  of  the  water  depth  and  the  cross¬ 
shore  and  alongshore  depth-averaged  and  near  bottom  velocities. 

The  developed  2D  and  quasi-3D  models  are  compared  with  the  laboratory 
data  for  obliquely  incident  regular  waves  of  Visser  (1991),  who  ensured  the  along¬ 
shore  uniformity  of  longshore  currents  using  a  pumping  system,  and  with  field  data 
of  Thornton  and  Guza  (1986)  for  that  the  assumption  of  alongshore  uniformity  ap¬ 
peared  reasonable.  To  assess  the  importance  of  the  dispersion  eflFects,  both  models 
are  compared  with  the  same  data  sets.  Comparisons  between  the  measured  and 
computed  cross-shore  variations  of  the  wave  height  and  setup  for  regular  waves  and 
the  root-mean-square  wave  height  for  irregular  waves  indicate  that  both  models  are 
capable  of  predicting  these  time-averaged  quantities.  The  dispersion  effects  on  the 
cross-shore  variations  of  the  wave  height  and  setup  are  shown  to  be  minor,  indicat¬ 
ing  that  the  cross-shore  dispersion  term  m  may  be  neglected  in  the  depth-integrated 
cross-shore  momentum  equation  as  anticipated  by  Kobayashi  and  Wurjanto  (1992). 

The  computed  oscillatory  components  of  the  cross-shore  and  alongshore  depth- 
averaged  velocities  become  the  maximum  near  the  still  water  shoreline  for  the  field 
data  of  Thornton  and  Guza  (1986),  contrary  to  the  conventional  assumption  of 
depth-limited  linear  breaking  waves  made  in  existing  time- averaged  models.  The 


causes  of  these  large  velocities  are  examined  using  the  time-averaged  wave  en¬ 
ergy  equation  and  filtering  the  computed  temporal  variations  in  the  high-  and 
low-frequency  bands.  Xhe  computed  energy  dissipation  rate  indicates  intense  wave 
breaking  immediately  seaward  of  the  still  water  shoreline  that  may  partly  explain 
the  computed  large  velocities  near  the  shoreline.  The  computed  high-  and  low- 
frequency  components  of  the  alongshore  velocity  variance  are  of  the  same  magni¬ 
tude  near  the  shoreline  even  in  the  absence  of  incident  low  frequency  waves  and  edge 
waves.  Kobayashi  et  al.  (1997)  derived  simple  relationships  between  the  free  surface 
and  cross-shore  velocity  statistics  using  a  linear  progressive  long-wave  theory.  The 
computed  values  of  the  normalized  free  surface  standard  deviation,  ^rms/^?  increase 
gradually  and  then  rapidly  near  the  shoreline  as  the  mean  water  depth  h  approaches 
zero.  The  computed  landward  increases  of  the  cross-shore  velocity  standard  devia¬ 
tion  and  the  undertow  magnitude  can  be  explained  using  these  relationships  with 
the  computed  cross-shore  variation  of  rj^jh. 

For  regular  waves,  the  2D  model  with  the  bottom  friction  factor  calibrated 
previously  for  swash  oscillations  predicts  the  magnitude  of  the  longshore  current 
but  cannot  reproduce  the  longshore  current  profile.  The  quasi-3D  model,  which 
includes  the  dispersion  effects,  improves  the  prediction  of  the  longshore  current 
profile  significantly.  For  irregular  waves,  both  models  using  the  bottom  friction 
factor  calibrated  by  Raubenheimer  et  al.  (1995),  predict  the  longshore  current 
profile  as  accurately  as  the  calibrated  time-averaged  model  by  Thornton  and  Guza 
(1986),  which  requires  much  less  computational  efforts.  The  advantage  of  the  time- 
dependent  models  can  not  be  confirmed  for  lack  of  swash  velocity  data.  The  3D 
model  computation  is  also  made  using  different  values  of  the  bottom  friction  factor 
to  improve  the  agreement  of  the  longshore  current.  The  bottom  friction  factor  is 
important  in  determining  the  magnitude  of  the  longshore  current  but  modifies  its 
profile  little.  In  summary,  the  dispersion  effects  on  the  longshore  current  profile  are 
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significant  for  regular  waves  but  secondary  for  irregular  waves,  especially  in  view  of 
the  uncertainties  associated  with  the  bottom  friction  factor.  Furthermore,  the  3D 
model  is  shown  to  predict  the  vertical  variations  of  the  longshore  current  inside  the 
surf  zone  but  not  outside  the  surf  zone  since  the  assumed  cubic  velocity  profile  may 
not  be  appropriate  outside  the  surf  zone. 

The  computed  results  using  the  2D  model  are  used  to  explain  the  very  low- 
frequency  component  of  the  computed  depth-averaged  alongshore  velocities  V  gen¬ 
erated  by  irregular  breaking  waves.  The  oscillatory  components  of  the  computed 
alongshore  velocity  are  analyzed  using  the  nonlinear  and  linear  alongshore  momen¬ 
tum  equations  together  with  the  computed  alongshore  gradient  of  the  free  surface 
elevation.  The  nonlinear  terms  associated  with  the  alongshore  momentum  flux  and 
bottom  shear  stress  are  shown  to  be  negligible  outside  the  surf  zone  and  reduce 
the  low-frequency  components  of  the  alongshore  velocity  in  the  breaker  and  inner 
surf  zones.  On  the  other  hand,  the  alongshore  free  surface  gradient  is  found  to  be 
negligible  in  driving  the  alongshore  oscillatory  fluid  motion  in  the  swash  zone  except 
for  the  very  low-frequency  range.  The  very  slow  oscillations  of  the  computed  along¬ 
shore  velocity  in  the  vicinity  of  the  breaker  zone  are  inferred  to  be  caused  by  the 
very  low-frequency  component  of  the  alongshore  free  surface  gradient  whose  origin 
is  uncertain  and  can  be  numerical. 

Finally,  the  3D  model  is  also  compared  with  the  DELILAH  field  data  for  a 
barred  beach  (Smith  et  al.  1993).  The  frequency  and  directional  spectra  were 
measured  at  the  8  m  depth.  The  seaward  boundary  of  the  numerical  model,  based 
on  the  assumption  of  shallow  water  waves,  is  taken  at  the  3.64  m  depth.  The 
shoaling  and  refraction  of  linear  directional  random  waves  is  computed  to  estimate 
the  frequency  spectrum  and  dominant  wave  direction  at  the  seaward  boundary. 
The  3D  model  underpredicts  the  root-mean-square  wave  height  somewhat  in  the 
bar  trough  region  and  predicts  a  peak  in  the  longshore  current  at  the  seaward 
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edge  of  the  bar  crest.  Under  the  assumption  of  alongshore  uniformity,  the  model 
cannot  explain  the  observed  broad  peak  in  the  longshore  current  in  the  bar  trough 
region.  Small  alongshore  variations  of  the  incident  wave  intensity  and  resulting  wave 
setup  are  shown  to  modify  the  longshore  current  profile  in  the  bar  trough  region 
significantly.  The  cross-shore  gradient  of  the  alongshore  radiation  stress  driving  the 
longshore  current  is  very  small  in  this  bar  trough  region.  The  alongshore  gradient 
of  wave  setup  is  shown  to  alter  the  force  driving  the  longshore  current  significantly 
in  this  region  and  produces  a  broad  peak  in  the  longshore  current.  Contrary  to 
the  computed  results  for  the  barred  beach,  the  longshore  current  profile  on  planar 
beaches  is  found  to  be  insensitive  to  the  alongshore  variations  of  the  incident  wave 
intensity  and  resulting  wave  setup.  As  a  result,  the  prediction  of  the  longshore 
current  profiles  on  barred  beaches  will  require  the  knowledge  of  small  alongshore 
variability  that  is  very  difficult  to  measure  accurately. 
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